diff --git a/DataFormats/L1THGCal/interface/HGCalCluster.h b/DataFormats/L1THGCal/interface/HGCalCluster.h index 99359f5f22bbf..fb84f17194102 100644 --- a/DataFormats/L1THGCal/interface/HGCalCluster.h +++ b/DataFormats/L1THGCal/interface/HGCalCluster.h @@ -8,7 +8,7 @@ namespace l1t { - + class HGCalCluster : public HGCalClusterT { public: diff --git a/DataFormats/L1THGCal/interface/HGCalClusterT.h b/DataFormats/L1THGCal/interface/HGCalClusterT.h index e0e07027cf9fd..1db937c448546 100644 --- a/DataFormats/L1THGCal/interface/HGCalClusterT.h +++ b/DataFormats/L1THGCal/interface/HGCalClusterT.h @@ -1,14 +1,15 @@ #ifndef DataFormats_L1Trigger_HGCalClusterT_h #define DataFormats_L1Trigger_HGCalClusterT_h +/* CMSSW */ #include "DataFormats/Common/interface/Ptr.h" -#include "DataFormats/Common/interface/PtrVector.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include "DataFormats/L1Trigger/interface/L1Candidate.h" #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h" #include "DataFormats/L1THGCal/interface/ClusterShapes.h" -#include "Math/Vector3D.h" +/* ROOT */ +#include "Math/Vector3D.h" namespace l1t { @@ -16,7 +17,7 @@ namespace l1t { public: - typedef typename edm::PtrVector::const_iterator const_iterator; + typedef typename std::vector>::const_iterator const_iterator; public: HGCalClusterT(){} @@ -43,63 +44,129 @@ namespace l1t { addConstituent(c); } - + ~HGCalClusterT() {}; - - const edm::PtrVector& constituents() const {return constituents_;} - const_iterator constituents_begin() const {return constituents_.begin();} - const_iterator constituents_end() const {return constituents_.end();} + + const std::vector>& constituents() const { return constituents_; } + const_iterator constituents_begin() const { return constituents_.begin(); } + const_iterator constituents_end() const { return constituents_.end(); } unsigned size() const { return constituents_.size(); } - void addConstituent( const edm::Ptr& c ) + void addConstituent( const edm::Ptr& c, bool updateCentre=true, float fraction=1. ) { - if( constituents_.empty() ) - { - detId_ = HGCalDetId(c->detId()); - seedMipPt_ = c->mipPt(); - } - /* update cluster positions */ - Basic3DVector constituentCentre( c->position() ); - Basic3DVector clusterCentre( centre_ ); + double cMipt = c->mipPt()*fraction; - clusterCentre = clusterCentre*mipPt_ + constituentCentre*c->mipPt(); - if( mipPt_ + c->mipPt()!=0 ) + if( constituents_.empty() ) { - clusterCentre /= ( mipPt_ + c->mipPt() ) ; + detId_ = HGCalDetId( c->detId() ); + seedMipPt_ = cMipt; + /* if the centre will not be dynamically calculated + the seed centre is considere as cluster centre */ + if( !updateCentre ) + { + centre_ = c->position(); + } } - centre_ = GlobalPoint( clusterCentre ); - if( clusterCentre.z()!=0 ) - { - centreProj_= GlobalPoint( clusterCentre / clusterCentre.z() ); + /* update cluster positions (IF requested) */ + if( updateCentre ){ + Basic3DVector constituentCentre( c->position() ); + Basic3DVector clusterCentre( centre_ ); + + clusterCentre = clusterCentre*mipPt_ + constituentCentre*cMipt; + if( (mipPt_ + cMipt ) > 0 ) + { + clusterCentre /= ( mipPt_ + cMipt ); + } + centre_ = GlobalPoint( clusterCentre ); + + if( clusterCentre.z()!=0 ) + { + centreProj_= GlobalPoint( clusterCentre / clusterCentre.z() ); + } } + /* update cluster energies */ - mipPt_ += c->mipPt(); + mipPt_ += cMipt; - int updatedPt = hwPt() + c->hwPt(); - setHwPt(updatedPt); + int updatedPt = hwPt() + (int)(c->hwPt()*fraction); + setHwPt( updatedPt ); math::PtEtaPhiMLorentzVector updatedP4 ( p4() ); - updatedP4 += c->p4(); + updatedP4 += (c->p4()*fraction); setP4( updatedP4 ); constituents_.push_back( c ); + constituentsFraction_.push_back( fraction ); } - - bool valid() const { return valid_;} - void setValid(bool valid) { valid_ = valid;} - + + void removeConstituent( const edm::Ptr& c, bool updateCentre=true ){ + + /* remove the pointer to c from the std::vector */ + double fraction=0; + bool constituentRemoved=false; + for( unsigned i=0; imipPt()*fraction; + if( updateCentre ){ + Basic3DVector constituentCentre( c->position() ); + Basic3DVector clusterCentre( centre_ ); + + clusterCentre = clusterCentre*mipPt_ - constituentCentre*cMipt; + if( (mipPt_ - cMipt ) > 0 ) + { + clusterCentre /= ( mipPt_ - cMipt ) ; + } + centre_ = GlobalPoint( clusterCentre ); + + if( clusterCentre.z() != 0 ) + { + centreProj_= GlobalPoint( clusterCentre / clusterCentre.z() ); + } + + } + + /* update cluster energies */ + mipPt_ -= cMipt; + + int updatedPt = hwPt() - ( c->hwPt()*fraction ); + setHwPt( updatedPt ); + + math::PtEtaPhiMLorentzVector updatedP4 ( p4() ); + updatedP4 -= ( c->p4()*fraction ); + setP4( updatedP4 ); + + } + + } + + bool valid() const { return valid_; } + void setValid(bool valid) { valid_ = valid; } + double mipPt() const { return mipPt_; } double seedMipPt() const { return seedMipPt_; } uint32_t detId() const { return detId_.rawId(); } + /* distance in 'cm' */ - double distance( const l1t::HGCalTriggerCell &tc ) const - { - return ( tc.position() - centre_ ).mag(); - } + double distance( const l1t::HGCalTriggerCell &tc ) const { return ( tc.position() - centre_ ).mag(); } const GlobalPoint& position() const { return centre_; } const GlobalPoint& centre() const { return centre_; } @@ -139,7 +206,37 @@ namespace l1t uint32_t subdetId() const {return detId_.subdetId();} uint32_t layer() const {return detId_.layer();} int32_t zside() const {return detId_.zside();} - + + + //shower shape + + int showerLength() const { return showerLength_; } + int coreShowerLength() const { return coreShowerLength_; } + int firstLayer() const { return firstLayer_; } + int maxLayer() const { return maxLayer_; } + float eMax() const { return eMax_; } + float sigmaEtaEtaMax() const { return sigmaEtaEtaMax_; } + float sigmaPhiPhiMax() const { return sigmaPhiPhiMax_; } + float sigmaEtaEtaTot() const { return sigmaEtaEtaTot_; } + float sigmaPhiPhiTot() const { return sigmaPhiPhiTot_; } + float sigmaZZ() const { return sigmaZZ_; } + float sigmaRRTot() const { return sigmaRRTot_; } + float sigmaRRMax() const { return sigmaRRMax_; } + float sigmaRRMean() const { return sigmaRRMean_; } + + void set_showerLength(int showerLength) { showerLength_ = showerLength;} + void set_coreShowerLength(int coreShowerLength) { coreShowerLength_ = coreShowerLength;} + void set_firstLayer(int firstLayer) { firstLayer_ = firstLayer;} + void set_maxLayer(int maxLayer) { maxLayer_ = maxLayer;} + void set_eMax(float eMax) { eMax_ = eMax;} + void set_sigmaEtaEtaMax(float sigmaEtaEtaMax) { sigmaEtaEtaMax_ = sigmaEtaEtaMax;} + void set_sigmaEtaEtaTot(float sigmaEtaEtaTot) { sigmaEtaEtaTot_ = sigmaEtaEtaTot;} + void set_sigmaPhiPhiMax(float sigmaPhiPhiMax) { sigmaPhiPhiMax_ = sigmaPhiPhiMax;} + void set_sigmaPhiPhiTot(float sigmaPhiPhiTot) { sigmaPhiPhiTot_ = sigmaPhiPhiTot;} + void set_sigmaRRMax(float sigmaRRMax) { sigmaRRMax_ = sigmaRRMax;} + void set_sigmaRRTot(float sigmaRRTot) { sigmaRRTot_ = sigmaRRTot;} + void set_sigmaRRMean(float sigmaRRMean) { sigmaRRMean_ = sigmaRRMean;} + void set_sigmaZZ(float sigmaZZ) { sigmaZZ_ = sigmaZZ;} /* operators */ bool operator<(const HGCalClusterT& cl) const {return mipPt() < cl.mipPt();} @@ -147,17 +244,37 @@ namespace l1t bool operator<=(const HGCalClusterT& cl) const { return !(cl>*this); } bool operator>=(const HGCalClusterT& cl) const { return !(cl<*this); } + private: - + bool valid_; - HGCalDetId detId_; - edm::PtrVector constituents_; + HGCalDetId detId_; + + std::vector> constituents_; /* ???? possibly change this in something like */ + std::vector constituentsFraction_; /* vector,float>> ???? */ + GlobalPoint centre_; GlobalPoint centreProj_; // centre projected onto the first HGCal layer double mipPt_; double seedMipPt_; + //shower shape + + int showerLength_; + int coreShowerLength_; + int firstLayer_; + int maxLayer_; + float eMax_; + float sigmaEtaEtaMax_; + float sigmaPhiPhiMax_; + float sigmaRRMax_; + float sigmaEtaEtaTot_; + float sigmaPhiPhiTot_; + float sigmaRRTot_; + float sigmaRRMean_; + float sigmaZZ_; + ClusterShapes shapes_; }; diff --git a/DataFormats/L1THGCal/src/classes_def.xml b/DataFormats/L1THGCal/src/classes_def.xml index ee8b31cc5815b..3ebfc6de72fd1 100644 --- a/DataFormats/L1THGCal/src/classes_def.xml +++ b/DataFormats/L1THGCal/src/classes_def.xml @@ -1,4 +1,3 @@ - @@ -21,16 +20,24 @@ - - + + + + + + - - + + + + + + diff --git a/L1Trigger/L1THGCal/data/panel_mapping_tdr_0.txt b/L1Trigger/L1THGCal/data/panel_mapping_tdr_0.txt new file mode 100644 index 0000000000000..5ddcff6c81229 --- /dev/null +++ b/L1Trigger/L1THGCal/data/panel_mapping_tdr_0.txt @@ -0,0 +1,468 @@ +7 33 +8 11 +9 1 +10 171 +11 161 +12 139 +13 129 +14 107 +15 97 +16 75 +17 65 +18 43 +19 33 +20 11 +21 4 +22 1 +23 171 +24 164 +25 161 +26 139 +27 132 +28 129 +29 107 +30 100 +31 97 +32 75 +33 68 +34 65 +35 43 +36 36 +37 33 +38 11 +39 4 +40 4 +41 1 +42 171 +43 164 +44 164 +45 161 +46 139 +47 132 +48 132 +49 129 +50 107 +51 100 +52 100 +53 97 +54 75 +55 68 +56 68 +57 65 +58 43 +59 36 +60 36 +61 34 +62 12 +63 14 +64 7 +65 5 +66 2 +67 172 +68 174 +69 167 +70 165 +71 162 +72 140 +73 142 +74 135 +75 133 +76 130 +77 108 +78 110 +79 103 +80 101 +81 98 +82 76 +83 78 +84 71 +85 69 +86 66 +87 44 +88 46 +89 39 +90 37 +91 34 +92 12 +93 14 +94 9 +95 7 +96 5 +97 2 +98 172 +99 174 +100 169 +101 167 +102 165 +103 162 +104 140 +105 142 +106 137 +107 135 +108 133 +109 130 +110 108 +111 110 +112 105 +113 103 +114 101 +115 98 +116 76 +117 78 +118 73 +119 71 +120 69 +121 66 +122 44 +123 46 +124 41 +125 39 +126 37 +127 34 +128 12 +129 14 +130 16 +131 9 +132 7 +133 5 +134 2 +135 172 +136 174 +137 176 +138 169 +139 167 +140 165 +141 162 +142 140 +143 142 +144 144 +145 137 +146 135 +147 133 +148 130 +149 108 +150 110 +151 112 +152 105 +153 103 +154 101 +155 98 +156 76 +157 78 +158 80 +159 73 +160 71 +161 69 +162 66 +163 44 +164 46 +165 48 +166 41 +167 39 +168 37 +169 35 +170 13 +171 15 +172 16 +173 17 +174 9 +175 8 +176 6 +177 3 +178 173 +179 175 +180 176 +181 177 +182 169 +183 168 +184 166 +185 163 +186 141 +187 143 +188 144 +189 145 +190 137 +191 136 +192 134 +193 131 +194 109 +195 111 +196 112 +197 113 +198 105 +199 104 +200 102 +201 99 +202 77 +203 79 +204 80 +205 81 +206 73 +207 72 +208 70 +209 67 +210 45 +211 47 +212 48 +213 49 +214 41 +215 40 +216 38 +217 35 +218 13 +219 15 +220 16 +221 17 +222 18 +223 10 +224 8 +225 6 +226 3 +227 173 +228 175 +229 176 +230 177 +231 178 +232 170 +233 168 +234 166 +235 163 +236 141 +237 143 +238 144 +239 145 +240 146 +241 138 +242 136 +243 134 +244 131 +245 109 +246 111 +247 112 +248 113 +249 114 +250 106 +251 104 +252 102 +253 99 +254 77 +255 79 +256 80 +257 81 +258 82 +259 74 +260 72 +261 70 +262 67 +263 45 +264 47 +265 48 +266 49 +267 50 +268 42 +269 40 +270 38 +271 35 +272 13 +273 15 +274 16 +275 17 +276 19 +277 18 +278 10 +279 8 +280 6 +281 3 +282 173 +283 175 +284 176 +285 177 +286 179 +287 178 +288 170 +289 168 +290 166 +291 163 +292 141 +293 143 +294 144 +295 145 +296 147 +297 146 +298 138 +299 136 +300 134 +301 131 +302 109 +303 111 +304 112 +305 113 +306 115 +307 114 +308 106 +309 104 +310 102 +311 99 +312 77 +313 79 +314 80 +315 81 +316 83 +317 82 +318 74 +319 72 +320 70 +321 67 +322 45 +323 47 +324 48 +325 49 +326 51 +327 50 +328 42 +329 40 +330 38 +331 35 +332 13 +333 15 +334 16 +335 17 +336 19 +337 19 +338 18 +339 10 +340 8 +341 6 +342 3 +343 173 +344 175 +345 176 +346 177 +347 179 +348 179 +349 178 +350 170 +351 168 +352 166 +353 163 +354 141 +355 143 +356 144 +357 145 +358 147 +359 147 +360 146 +361 138 +362 136 +363 134 +364 131 +365 109 +366 111 +367 112 +368 113 +369 115 +370 115 +371 114 +372 106 +373 104 +374 102 +375 99 +376 77 +377 79 +378 80 +379 81 +380 83 +381 83 +382 82 +383 74 +384 72 +385 70 +386 67 +387 45 +388 47 +389 48 +390 49 +391 51 +392 51 +393 50 +394 42 +395 40 +396 38 +399 15 +400 16 +401 17 +402 19 +403 19 +404 19 +405 18 +406 10 +407 8 +411 175 +412 176 +413 177 +414 179 +415 179 +416 179 +417 178 +418 170 +419 168 +423 143 +424 144 +425 145 +426 147 +427 147 +428 147 +429 146 +430 138 +431 136 +435 111 +436 112 +437 113 +438 115 +439 115 +440 115 +441 114 +442 106 +443 104 +447 79 +448 80 +449 81 +450 83 +451 83 +452 83 +453 82 +454 74 +455 72 +459 47 +460 48 +461 49 +462 51 +463 51 +464 51 +465 50 +466 42 +467 40 +474 20 +475 20 +476 20 +477 20 +487 180 +488 180 +489 180 +490 180 +500 148 +501 148 +502 148 +503 148 +513 116 +514 116 +515 116 +516 116 +526 84 +527 84 +528 84 +529 84 +539 52 +540 52 +541 52 +542 52 diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h b/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h index 13f7902926ed9..69f7b213b79e0 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h @@ -62,6 +62,7 @@ class HGCalTriggerGeometryBase virtual bool validTriggerCell( const unsigned trigger_cell_id) const = 0; virtual bool disconnectedModule(const unsigned module_id) const = 0; + virtual unsigned triggerLayer(const unsigned id) const = 0; protected: void setCaloGeometry(const edm::ESHandle& geom) {calo_geometry_=geom;} diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryGenericMapping.h b/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryGenericMapping.h index cd1669036358e..9f19d497d64d7 100644 --- a/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryGenericMapping.h +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerGeometryGenericMapping.h @@ -139,6 +139,7 @@ class HGCalTriggerGeometryGenericMapping : public HGCalTriggerGeometryBase { virtual bool validTriggerCell( const unsigned trigger_cell_det_id ) const override final; virtual bool disconnectedModule(const unsigned module_id) const override final; + virtual unsigned triggerLayer(const unsigned id) const override final; protected: geom_map cells_to_trigger_cells_; diff --git a/L1Trigger/L1THGCal/interface/HGCalTriggerTools.h b/L1Trigger/L1THGCal/interface/HGCalTriggerTools.h new file mode 100644 index 0000000000000..0c4960d21fee2 --- /dev/null +++ b/L1Trigger/L1THGCal/interface/HGCalTriggerTools.h @@ -0,0 +1,71 @@ +#ifndef __L1Trigger_L1THGCal_HGCalTriggerTools_h__ +#define __L1Trigger_L1THGCal_HGCalTriggerTools_h__ + +/** \class HGCalTriggerTools + * Tools for handling HGCal trigger det-ID: in the current version + * of trhe HGCAL simulation only HGCalDetId for the TriggerCells (TC) + * are used and not HcalDetId as in the offline! + * As a consequence the class assumes that only DetIds of the first kind are used in the getTC* methods + * NOTE: this uses the trigger geometry hence would give wrong results + * when used for offline reco!!!! + * + * \author G. Cerminara (CERN), heavily "inspired" by HGCalRechHitTools ;) + */ + + +#include +#include +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h" + +class HGCalTriggerGeometryBase; +class DetId; + + + + + +namespace edm { + class Event; + class EventSetup; +} + + class HGCalTriggerTools { + public: + HGCalTriggerTools() : geom_(nullptr), + fhOffset_(0), + bhOffset_(0) {} + ~HGCalTriggerTools() {} + + void setEventSetup(const edm::EventSetup&); + GlobalPoint getTCPosition(const DetId& id) const; + unsigned int getLayerWithOffset(const DetId&) const; + // unsigned int getLayer(ForwardSubdetector type) const; + unsigned int getLayer(const DetId&) const; + + // 4-vector helper functions using GlobalPoint + float getEta(const GlobalPoint& position, const float& vertex_z = 0.) const; + float getPhi(const GlobalPoint& position) const; + float getPt(const GlobalPoint& position, const float& hitEnergy, const float& vertex_z = 0.) const; + + // 4-vector helper functions using DetId + float getTCEta(const DetId& id, const float& vertex_z = 0.) const; + float getTCPhi(const DetId& id) const; + float getTCPt(const DetId& id, const float& hitEnergy, const float& vertex_z = 0.) const; + + inline const HGCalTriggerGeometryBase * getTriggerGeometry() const {return geom_;}; + unsigned int lastLayerEE() const {return fhOffset_;} + unsigned int lastLayerFH() const {return bhOffset_;} + + float getLayerZ(const unsigned& layerWithOffset) const; + float getLayerZ(const int& subdet, const unsigned& layer) const; + + + + private: + const HGCalTriggerGeometryBase* geom_; + unsigned int fhOffset_, bhOffset_; + }; + + +#endif diff --git a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalClusteringImpl.h b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalClusteringImpl.h index 689970f767704..7d05a2aca3749 100644 --- a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalClusteringImpl.h +++ b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalClusteringImpl.h @@ -10,6 +10,10 @@ #include "DataFormats/L1THGCal/interface/HGCalCluster.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" + +bool distanceSorter (pair,float> i,pair,float> j) { return (i.second & triggerCellsPtrs, - l1t::HGCalClusterBxCollection & clusters + void clusterizeDR( const std::vector> & triggerCellsPtrs, + l1t::HGCalClusterBxCollection & clusters ); /* NN-algorithms */ @@ -44,11 +48,18 @@ class HGCalClusteringImpl{ const HGCalTriggerGeometryBase & triggerGeometry ); - void clusterizeNN( const edm::PtrVector & triggerCellsPtrs, + void clusterizeNN( const std::vector> & triggerCellsPtrs, l1t::HGCalClusterBxCollection & clusters, const HGCalTriggerGeometryBase & triggerGeometry ); + /* FW-algorithms */ + void clusterizeDRNN( const std::vector> & triggerCellsPtrs, + l1t::HGCalClusterBxCollection & clusters, + const HGCalTriggerGeometryBase & triggerGeometry + ); + + private: @@ -58,9 +69,20 @@ class HGCalClusteringImpl{ double scintillatorTriggerCellThreshold_; double dr_; std::string clusteringAlgorithmType_; - void triggerCellReshuffling( const edm::PtrVector & triggerCellsPtrs, + double calibSF_; + std::vector layerWeights_; + bool applyLayerWeights_; + + void triggerCellReshuffling( const std::vector> & triggerCellsPtrs, std::array>, kLayers_>, kNSides_> & reshuffledTriggerCells ); + bool areTCneighbour( uint32_t detIDa, uint32_t detIDb, const HGCalTriggerGeometryBase & triggerGeometry ); + + void removeUnconnectedTCinCluster( l1t::HGCalCluster & cluster, + const HGCalTriggerGeometryBase & triggerGeometry + ); + + void calibratePt( l1t::HGCalCluster & cluster ); }; diff --git a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalMulticlusteringImpl.h b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalMulticlusteringImpl.h index 4c2198f341d5a..8c3e25ffe1593 100644 --- a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalMulticlusteringImpl.h +++ b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalMulticlusteringImpl.h @@ -6,6 +6,10 @@ #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" +#include "L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h" + + class HGCalMulticlusteringImpl{ public: @@ -16,14 +20,35 @@ class HGCalMulticlusteringImpl{ const l1t::HGCalMulticluster & mclu, double dR ) const; - void clusterize( const edm::PtrVector & clustersPtr, - l1t::HGCalMulticlusterBxCollection & multiclusters); + void clusterizeDR( const std::vector> & clustersPtr, + l1t::HGCalMulticlusterBxCollection & multiclusters, + const HGCalTriggerGeometryBase & triggerGeometry + ); + + void clusterizeDBSCAN( const std::vector> & clustersPtr, + l1t::HGCalMulticlusterBxCollection & multiclusters, + const HGCalTriggerGeometryBase & triggerGeometry + ); private: + + void findNeighbor( const std::vector>& rankedList, + unsigned int searchInd, + const std::vector> & clustersPtr, + std::vector& neigbors); double dr_; double ptC3dThreshold_; - double calibSF_; + string multiclusterAlgoType_; + double distDbscan_ = 0.005; + unsigned minNDbscan_ = 3; + + HGCalShowerShape shape_; + + static const int kLayersEE_=28; + static const int kLayersFH_=12; + static const int kLayersBH_=12; + }; #endif diff --git a/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h new file mode 100644 index 0000000000000..2feb0b31ed03c --- /dev/null +++ b/L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h @@ -0,0 +1,58 @@ +#ifndef __L1Trigger_L1THGCal_HGCALSHOWERSHAPE_h__ +#define __L1Trigger_L1THGCal_HGCALSHOWERSHAPE_h__ +#include +#include +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" +#include "DataFormats/Math/interface/LorentzVector.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" + +class HGCalShowerShape{ + + public: + typedef math::XYZTLorentzVector LorentzVector; + + HGCalShowerShape(){} + + ~HGCalShowerShape(){} + + int firstLayer(const l1t::HGCalMulticluster& c3d) const; + int lastLayer(const l1t::HGCalMulticluster& c3d) const; + int maxLayer(const l1t::HGCalMulticluster& c3d) const; + int showerLength(const l1t::HGCalMulticluster& c3d) const {return lastLayer(c3d)-firstLayer(c3d)+1; }//in number of layers + // Maximum number of consecutive layers in the cluster + int coreShowerLength(const l1t::HGCalMulticluster& c3d, const HGCalTriggerGeometryBase& triggerGeometry) const; + + float eMax(const l1t::HGCalMulticluster& c3d) const; + + float sigmaZZ(const l1t::HGCalMulticluster& c3d) const; + + float sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const; + float sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const; + float sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const; + + float sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const; + float sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const; + float sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const; + + float sigmaRRTot(const l1t::HGCalMulticluster& c3d) const; + float sigmaRRTot(const l1t::HGCalCluster& c2d) const; + float sigmaRRMax(const l1t::HGCalMulticluster& c3d) const; + float sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius=5.) const; + + private: + + float meanX(const std::vector >& energy_X_tc) const; + float sigmaXX(const std::vector >& energy_X_tc, const float X_cluster) const; + float sigmaPhiPhi(const std::vector >& energy_phi_tc, const float phi_cluster) const; + + static const int kLayersEE_=28; + static const int kLayersFH_=12; + static const int kLayersBH_=12; + int HGC_layer(const uint32_t subdet, const uint32_t layer) const; + + +}; + + +#endif + diff --git a/L1Trigger/L1THGCal/plugins/BuildFile.xml b/L1Trigger/L1THGCal/plugins/BuildFile.xml index dc9a04cbc48d4..8782ce8325ae4 100644 --- a/L1Trigger/L1THGCal/plugins/BuildFile.xml +++ b/L1Trigger/L1THGCal/plugins/BuildFile.xml @@ -11,7 +11,7 @@ - + @@ -27,10 +27,18 @@ - + + + + + + + + + diff --git a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCClusterAlgo.cc b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCClusterAlgo.cc index 825cd76676de3..a738c0aa0f466 100644 --- a/L1Trigger/L1THGCal/plugins/be_algorithms/HGCClusterAlgo.cc +++ b/L1Trigger/L1THGCal/plugins/be_algorithms/HGCClusterAlgo.cc @@ -26,7 +26,12 @@ class HGCClusterAlgo : public Algorithm private: enum ClusterType{ dRC2d, - NNC2d + NNC2d, + dRNNC2d + }; + enum MulticlusterType{ + dRC3d, + DBSCANC3d }; public: @@ -38,20 +43,32 @@ class HGCClusterAlgo : public Algorithm multicluster_product_( new l1t::HGCalMulticlusterBxCollection ), calibration_( conf.getParameterSet("calib_parameters") ), clustering_( conf.getParameterSet("C2d_parameters") ), - multiclustering_( conf.getParameterSet("C3d_parameters" ) ) + multiclustering_( conf.getParameterSet("C3d_parameters" ) ), + triggercell_threshold_silicon_( conf.getParameter("triggercell_threshold_silicon") ), + triggercell_threshold_scintillator_( conf.getParameter("triggercell_threshold_scintillator") ) { - clustering_threshold_silicon_ = conf.getParameterSet("C2d_parameters").getParameter("clustering_threshold_silicon"); - clustering_threshold_scintillator_ = conf.getParameterSet("C2d_parameters").getParameter("clustering_threshold_scintillator"); std::string type(conf.getParameterSet("C2d_parameters").getParameter("clusterType")); if(type=="dRC2d"){ clusteringAlgorithmType_ = dRC2d; }else if(type=="NNC2d"){ clusteringAlgorithmType_ = NNC2d; + }else if(type=="dRNNC2d"){ + clusteringAlgorithmType_ = dRNNC2d; }else { edm::LogWarning("ParameterError") << "Unknown clustering type '" << type << "'. Using nearest neighbor NNC2d instead.\n"; clusteringAlgorithmType_ = NNC2d; } + std::string typeMulticluster(conf.getParameterSet("C3d_parameters").getParameter("type_multicluster")); + if(typeMulticluster=="dRC3d"){ + multiclusteringAlgoType_ = dRC3d; + }else if(typeMulticluster=="DBSCANC3d"){ + multiclusteringAlgoType_ = DBSCANC3d; + }else { + edm::LogWarning("ParameterError") << "Unknown Multiclustering type '" << typeMulticluster + << "'. Using Cone Algorithm instead.\n"; + multiclusteringAlgoType_ = dRC3d; + } } @@ -95,8 +112,9 @@ class HGCClusterAlgo : public Algorithm /* algorithm type */ ClusterType clusteringAlgorithmType_; - double clustering_threshold_silicon_; - double clustering_threshold_scintillator_; + double triggercell_threshold_silicon_; + double triggercell_threshold_scintillator_; + MulticlusterType multiclusteringAlgoType_; }; @@ -122,8 +140,8 @@ void HGCClusterAlgo::run(const l1t::HGCFETriggerDigiCollection & c { l1t::HGCalTriggerCell calibratedtriggercell( triggercell ); calibration_.calibrateInGeV( calibratedtriggercell); - double clustering_threshold = (triggercell.subdetId()==HGCHEB ? clustering_threshold_scintillator_ : clustering_threshold_silicon_); - if(calibratedtriggercell.mipPt()push_back( 0, calibratedtriggercell ); } @@ -140,7 +158,7 @@ void HGCClusterAlgo::run(const l1t::HGCFETriggerDigiCollection & c triggerCellsHandle = evt.put( std::move( trgcell_product_ ), "calibratedTriggerCells"); /* create a persistent vector of pointers to the trigger-cells */ - edm::PtrVector triggerCellsPtrs; + std::vector> triggerCellsPtrs; for( unsigned i = 0; i < triggerCellsHandle->size(); ++i ) { edm::Ptr ptr(triggerCellsHandle,i); triggerCellsPtrs.push_back(ptr); @@ -154,6 +172,9 @@ void HGCClusterAlgo::run(const l1t::HGCFETriggerDigiCollection & c case NNC2d: clustering_.clusterizeNN( triggerCellsPtrs, *cluster_product_, *triggerGeometry_ ); break; + case dRNNC2d: + clustering_.clusterizeDRNN( triggerCellsPtrs, *cluster_product_, *triggerGeometry_ ); + break; default: // Should not happen, clustering type checked in constructor break; @@ -163,18 +184,30 @@ void HGCClusterAlgo::run(const l1t::HGCFETriggerDigiCollection & c clustersHandle = evt.put( std::move( cluster_product_ ), "cluster2D"); /* create a persistent vector of pointers to the trigger-cells */ - edm::PtrVector clustersPtrs; + std::vector> clustersPtrs; for( unsigned i = 0; i < clustersHandle->size(); ++i ) { edm::Ptr ptr(clustersHandle,i); clustersPtrs.push_back(ptr); } - /* call to multiclustering */ - multiclustering_.clusterize( clustersPtrs, *multicluster_product_ ); + /* call to multiclustering and compute shower shape*/ + switch(multiclusteringAlgoType_){ + case dRC3d : + multiclustering_.clusterizeDR( clustersPtrs, *multicluster_product_, *triggerGeometry_); + break; + case DBSCANC3d: + multiclustering_.clusterizeDBSCAN( clustersPtrs, *multicluster_product_, *triggerGeometry_); + break; + default: + // Should not happen, clustering type checked in constructor + break; + } /* retrieve the orphan handle to the multiclusters collection and put the collection in the event */ multiclustersHandle = evt.put( std::move( multicluster_product_ ), "cluster3D"); + + } typedef HGCClusterAlgo HGCClusterAlgoBestChoice; diff --git a/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexImp2.cc b/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexImp2.cc index 42f18985e86f8..476683c8483e1 100644 --- a/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexImp2.cc +++ b/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexImp2.cc @@ -36,6 +36,7 @@ class HGCalTriggerGeometryHexImp2 : public HGCalTriggerGeometryBase virtual bool validTriggerCell( const unsigned ) const override final; virtual bool disconnectedModule(const unsigned) const override final; + virtual unsigned triggerLayer(const unsigned) const override final; private: edm::FileInPath l1tCellsMapping_; @@ -773,6 +774,14 @@ disconnectedModule(const unsigned module_id) const return false; } + +unsigned +HGCalTriggerGeometryHexImp2:: +triggerLayer(const unsigned id) const +{ + return HGCalDetId(id).layer(); +} + bool HGCalTriggerGeometryHexImp2:: validTriggerCellFromCells(const unsigned trigger_cell_id) const diff --git a/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexLayerBasedImp1.cc b/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexLayerBasedImp1.cc index f0e2fe0212095..23ec1375b7198 100644 --- a/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexLayerBasedImp1.cc +++ b/L1Trigger/L1THGCal/plugins/geometries/HGCalTriggerGeometryHexLayerBasedImp1.cc @@ -37,6 +37,7 @@ class HGCalTriggerGeometryHexLayerBasedImp1 : public HGCalTriggerGeometryBase virtual bool validTriggerCell( const unsigned ) const override final; virtual bool disconnectedModule(const unsigned) const override final; + virtual unsigned triggerLayer(const unsigned) const override final; private: edm::FileInPath l1tCellsMapping_; @@ -65,8 +66,15 @@ class HGCalTriggerGeometryHexLayerBasedImp1 : public HGCalTriggerGeometryBase std::unordered_map>> trigger_cell_neighbors_; std::unordered_map>> trigger_cell_neighbors_bh_; - // Disconnected modules + // Disconnected modules and layers std::unordered_set disconnected_modules_; + std::unordered_set disconnected_layers_; + std::vector trigger_layers_; + + // layer offsets + unsigned fhOffset_; + unsigned bhOffset_; + unsigned totalLayers_; void fillMaps(); void fillNeighborMaps(const edm::FileInPath&, std::unordered_map>>&); @@ -81,6 +89,8 @@ class HGCalTriggerGeometryHexLayerBasedImp1 : public HGCalTriggerGeometryBase unsigned packIetaIphi(unsigned ieta, unsigned iphi) const; void unpackWaferCellId(unsigned wafer_cell, unsigned& wafer, unsigned& cell) const; void unpackIetaIphi(unsigned ieta_iphi, unsigned& ieta, unsigned& iphi) const; + + unsigned layerWithOffset(unsigned) const; }; @@ -95,6 +105,8 @@ HGCalTriggerGeometryHexLayerBasedImp1(const edm::ParameterSet& conf): { std::vector tmp_vector = conf.getParameter>("DisconnectedModules"); std::move(tmp_vector.begin(), tmp_vector.end(), std::inserter(disconnected_modules_, disconnected_modules_.end())); + tmp_vector = conf.getParameter>("DisconnectedLayers"); + std::move(tmp_vector.begin(), tmp_vector.end(), std::inserter(disconnected_layers_, disconnected_layers_.end())); } void @@ -118,6 +130,24 @@ HGCalTriggerGeometryHexLayerBasedImp1:: initialize(const edm::ESHandle& calo_geometry) { setCaloGeometry(calo_geometry); + fhOffset_ = eeTopology().dddConstants().layers(true); + bhOffset_ = fhOffset_ + fhTopology().dddConstants().layers(true); + totalLayers_ = bhOffset_ + bhTopology().dddConstants()->getMaxDepth(1); + trigger_layers_.resize(totalLayers_+1); + unsigned trigger_layer = 0; + for(unsigned layer=0; layer=trigger_layers_.size()) return 0; + return trigger_layers_[layer]; } bool @@ -761,6 +803,28 @@ detIdWaferType(unsigned subdet, short wafer) const } +unsigned +HGCalTriggerGeometryHexLayerBasedImp1:: +layerWithOffset(unsigned id) const +{ + HGCalDetId detid(id); + unsigned layer = 0; + switch(detid.subdetId()) + { + case ForwardSubdetector::HGCEE: + layer = detid.layer(); + break; + case ForwardSubdetector::HGCHEF: + layer = fhOffset_ + detid.layer(); + break; + case ForwardSubdetector::HGCHEB: + layer = bhOffset_ + detid.layer(); + break; + }; + return layer; +} + + DEFINE_EDM_PLUGIN(HGCalTriggerGeometryFactory, HGCalTriggerGeometryHexLayerBasedImp1, diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGen.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGen.cc index e229c617adf46..b17ff2a045fb2 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGen.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGen.cc @@ -1,6 +1,113 @@ #include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" +#include "DataFormats/GeometrySurface/interface/Plane.h" +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "DataFormats/HcalDetId/interface/HcalDetId.h" + + #include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h" +#include "TrackPropagation/RungeKutta/interface/RKPropagatorInS.h" +#include "FastSimulation/Event/interface/FSimEvent.h" +#include "FastSimulation/Particle/interface/ParticleTable.h" +#include "FastSimulation/CaloGeometryTools/interface/Transform3DPJ.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" + +// NOTE: most of this code is borrowed by https://github.com/CMS-HGCAL/reco-ntuples +// kudos goes to the original authors. Ideally the 2 repos should be merged since they share part of the use case +#include +#include + +namespace HGCal_helpers { + + class Coordinates { + public: + Coordinates() : x(0), y(0), z(0), eta(0), phi(0) {} + float x, y, z, eta, phi; + inline math::XYZTLorentzVectorD toVector() { return math::XYZTLorentzVectorD(x, y, z, 0); } + }; + + + class SimpleTrackPropagator { + public: + SimpleTrackPropagator(const MagneticField *f) + : field_(f), prod_(field_, alongMomentum), absz_target_(0) { + ROOT::Math::SMatrixIdentity id; + AlgebraicSymMatrix55 C(id); + const float uncert = 0.001; + C *= uncert; + err_ = CurvilinearTrajectoryError(C); + } + void setPropagationTargetZ(const float &z); + + bool propagate(const double px, const double py, const double pz, const double x, const double y, + const double z, const float charge, Coordinates &coords) const; + + bool propagate(const math::XYZTLorentzVectorD &momentum, const math::XYZTLorentzVectorD &position, + const float charge, Coordinates &coords) const; + + private: + SimpleTrackPropagator() : field_(0), prod_(field_, alongMomentum), absz_target_(0) {} + const RKPropagatorInS &RKProp() const { return prod_.propagator; } + Plane::PlanePointer targetPlaneForward_, targetPlaneBackward_; + const MagneticField *field_; + CurvilinearTrajectoryError err_; + defaultRKPropagator::Product prod_; + float absz_target_; + }; + void SimpleTrackPropagator::setPropagationTargetZ(const float &z) { + targetPlaneForward_ = Plane::build(Plane::PositionType(0, 0, std::abs(z)), Plane::RotationType()); + targetPlaneBackward_ = + Plane::build(Plane::PositionType(0, 0, -std::abs(z)), Plane::RotationType()); + absz_target_ = std::abs(z); + } + bool SimpleTrackPropagator::propagate(const double px, const double py, const double pz, + const double x, const double y, const double z, + const float charge, Coordinates &output) const { + output = Coordinates(); + + typedef TrajectoryStateOnSurface TSOS; + GlobalPoint startingPosition(x, y, z); + GlobalVector startingMomentum(px, py, pz); + Plane::PlanePointer startingPlane = + Plane::build(Plane::PositionType(x, y, z), Plane::RotationType()); + TSOS startingStateP( + GlobalTrajectoryParameters(startingPosition, startingMomentum, charge, field_), err_, + *startingPlane); + + TSOS trackStateP; + if (pz > 0) { + trackStateP = RKProp().propagate(startingStateP, *targetPlaneForward_); + } else { + trackStateP = RKProp().propagate(startingStateP, *targetPlaneBackward_); + } + if (trackStateP.isValid()) { + output.x = trackStateP.globalPosition().x(); + output.y = trackStateP.globalPosition().y(); + output.z = trackStateP.globalPosition().z(); + output.phi = trackStateP.globalPosition().phi(); + output.eta = trackStateP.globalPosition().eta(); + return true; + } + return false; + } + + bool SimpleTrackPropagator::propagate(const math::XYZTLorentzVectorD &momentum, + const math::XYZTLorentzVectorD &position, const float charge, + Coordinates &output) const { + return propagate(momentum.px(), momentum.py(), momentum.pz(), position.x(), position.y(), + position.z(), charge, output); + } + +} // HGCal_helpers class HGCalTriggerNtupleGen : public HGCalTriggerNtupleBase @@ -12,17 +119,87 @@ class HGCalTriggerNtupleGen : public HGCalTriggerNtupleBase virtual void initialize(TTree&, const edm::ParameterSet&, edm::ConsumesCollector&&) override final; virtual void fill(const edm::Event&, const edm::EventSetup& ) override final; + typedef ROOT::Math::Transform3DPJ Transform3D; + typedef ROOT::Math::Transform3DPJ::Point Point; + + enum ReachHGCal { + notReach = 0, + outsideEESurface = 1, + onEESurface = 2 + }; + private: virtual void clear() override final; edm::EDGetToken gen_token_; + edm::EDGetToken gen_PU_token_; + int gen_n_; - std::vector gen_id_; - std::vector gen_status_; - std::vector gen_energy_; - std::vector gen_pt_; + int gen_PUNumInt_; + float gen_TrueNumInt_; + + float vtx_x_; + float vtx_y_; + float vtx_z_; + + + //////////////////// + // GenParticles + // + std::vector genpart_eta_; + std::vector genpart_phi_; + std::vector genpart_pt_; + std::vector genpart_energy_; + std::vector genpart_dvx_; + std::vector genpart_dvy_; + std::vector genpart_dvz_; + std::vector genpart_ovx_; + std::vector genpart_ovy_; + std::vector genpart_ovz_; + std::vector genpart_exx_; + std::vector genpart_exy_; + std::vector genpart_mother_; + std::vector genpart_exphi_; + std::vector genpart_exeta_; + std::vector genpart_fbrem_; + std::vector genpart_pid_; + std::vector genpart_gen_; + std::vector genpart_reachedEE_; + std::vector genpart_fromBeamPipe_; + std::vector> genpart_posx_; + std::vector> genpart_posy_; + std::vector> genpart_posz_; + + //////////////////// + // reco::GenParticles + // std::vector gen_eta_; std::vector gen_phi_; + std::vector gen_pt_; + std::vector gen_energy_; + std::vector gen_charge_; + std::vector gen_pdgid_; + std::vector gen_status_; + std::vector> gen_daughters_; + + + // -------convenient tool to deal with simulated tracks + std::unique_ptr mySimEvent_; + + //std::vector dEdXWeights_; + //std::vector invThicknessCorrection_; + + // and also the magnetic field + const MagneticField *aField_; + + HGCalTriggerTools triggerTools_; + + // edm::EDGetTokenT > genParticles_; + edm::EDGetToken simTracks_token_; + edm::EDGetToken simVertices_token_; + edm::EDGetToken hepmcev_token_; + + }; DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, @@ -35,46 +212,239 @@ HGCalTriggerNtupleGen(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf { } + void HGCalTriggerNtupleGen:: initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector) { + edm::ParameterSet particleFilter_(conf.getParameter("particleFilter")); + mySimEvent_ = std::make_unique(particleFilter_); + gen_token_ = collector.consumes(conf.getParameter("GenParticles")); + gen_PU_token_ = collector.consumes>(conf.getParameter("GenPU")); tree.Branch("gen_n", &gen_n_, "gen_n/I"); - tree.Branch("gen_id", &gen_id_); - tree.Branch("gen_status", &gen_status_); - tree.Branch("gen_energy", &gen_energy_); - tree.Branch("gen_pt", &gen_pt_); + tree.Branch("gen_PUNumInt", &gen_PUNumInt_ ,"gen_PUNumInt/I"); + tree.Branch("gen_TrueNumInt", &gen_TrueNumInt_ ,"gen_TrueNumInt/F"); + + hepmcev_token_ = collector.consumes(edm::InputTag("generatorSmeared")); + + simTracks_token_ = collector.consumes>(edm::InputTag("g4SimHits")); + simVertices_token_ = collector.consumes>(edm::InputTag("g4SimHits")); + + tree.Branch("vtx_x", &vtx_x_); + tree.Branch("vtx_y", &vtx_y_); + tree.Branch("vtx_z", &vtx_z_); + + tree.Branch("gen_eta", &gen_eta_); tree.Branch("gen_phi", &gen_phi_); + tree.Branch("gen_pt", &gen_pt_); + tree.Branch("gen_energy", &gen_energy_); + tree.Branch("gen_charge", &gen_charge_); + tree.Branch("gen_pdgid", &gen_pdgid_); + tree.Branch("gen_status", &gen_status_); + tree.Branch("gen_daughters", &gen_daughters_); + + tree.Branch("genpart_eta", &genpart_eta_); + tree.Branch("genpart_phi", &genpart_phi_); + tree.Branch("genpart_pt", &genpart_pt_); + tree.Branch("genpart_energy", &genpart_energy_); + tree.Branch("genpart_dvx", &genpart_dvx_); + tree.Branch("genpart_dvy", &genpart_dvy_); + tree.Branch("genpart_dvz", &genpart_dvz_); + tree.Branch("genpart_ovx", &genpart_ovx_); + tree.Branch("genpart_ovy", &genpart_ovy_); + tree.Branch("genpart_ovz", &genpart_ovz_); + tree.Branch("genpart_mother", &genpart_mother_); + tree.Branch("genpart_exphi", &genpart_exphi_); + tree.Branch("genpart_exeta", &genpart_exeta_); + tree.Branch("genpart_exx", &genpart_exx_); + tree.Branch("genpart_exy", &genpart_exy_); + tree.Branch("genpart_fbrem", &genpart_fbrem_); + tree.Branch("genpart_pid", &genpart_pid_); + tree.Branch("genpart_gen", &genpart_gen_); + tree.Branch("genpart_reachedEE", &genpart_reachedEE_); + tree.Branch("genpart_fromBeamPipe", &genpart_fromBeamPipe_); + tree.Branch("genpart_posx", &genpart_posx_); + tree.Branch("genpart_posy", &genpart_posy_); + tree.Branch("genpart_posz", &genpart_posz_); + } void HGCalTriggerNtupleGen:: -fill(const edm::Event& e, const edm::EventSetup& es) +fill(const edm::Event& iEvent, const edm::EventSetup& es) { - edm::Handle gen_particles_h; - e.getByToken(gen_token_, gen_particles_h); - const reco::GenParticleCollection& gen_particles = *gen_particles_h; - clear(); - gen_n_ = gen_particles.size(); - gen_id_.reserve(gen_n_); - gen_status_.reserve(gen_n_); - gen_energy_.reserve(gen_n_); - gen_pt_.reserve(gen_n_); - gen_eta_.reserve(gen_n_); - gen_phi_.reserve(gen_n_); - for(const auto& particle : gen_particles) + + edm::Handle > PupInfo_h; + iEvent.getByToken(gen_PU_token_, PupInfo_h); + const std::vector< PileupSummaryInfo >& PupInfo = *PupInfo_h; + + + // FIXME: this part could go in begin run + edm::ESHandle pdt; + es.getData(pdt); + mySimEvent_->initializePdt(&(*pdt)); + + triggerTools_.setEventSetup(es); + + edm::ESHandle magfield; + es.get().get(magfield); + aField_ = &(*magfield); + // up to here...could go in the beginRun + + // This balck magic is needed to use the mySimEvent_ + ParticleTable::Sentry ptable(mySimEvent_->theTable()); + edm::Handle hevH; + edm::Handle> simTracksHandle; + edm::Handle> simVerticesHandle; + + iEvent.getByToken(hepmcev_token_, hevH); + iEvent.getByToken(simTracks_token_, simTracksHandle); + iEvent.getByToken(simVertices_token_, simVerticesHandle); + mySimEvent_->fill(*simTracksHandle, *simVerticesHandle); + + HepMC::GenVertex *primaryVertex = *(hevH)->GetEvent()->vertices_begin(); + const float mm2cm = 0.1; + vtx_x_ = primaryVertex->position().x() * mm2cm; // to put in official units + vtx_y_ = primaryVertex->position().y() * mm2cm; + vtx_z_ = primaryVertex->position().z() * mm2cm; + Point sim_pv(vtx_x_, vtx_y_, vtx_z_); + + + HGCal_helpers::SimpleTrackPropagator toHGCalPropagator(aField_); + toHGCalPropagator.setPropagationTargetZ(triggerTools_.getLayerZ(1)); + std::vector allselectedgentracks; + const float eeInnerRadius = 25.; + const float eeOuterRadius = 160.; + unsigned int npart = mySimEvent_->nTracks(); + for (unsigned int i = 0; i < npart; ++i) { + std::vector xp, yp, zp; + FSimTrack &myTrack(mySimEvent_->track(i)); + math::XYZTLorentzVectorD vtx(0, 0, 0, 0); + + int reachedEE = ReachHGCal::notReach; // compute the extrapolations for the particles reaching EE + // and for the gen particles + double fbrem = -1; + + if (std::abs(myTrack.vertex().position().z()) >= triggerTools_.getLayerZ(1)) continue; + + const unsigned nlayers = 52; + if (myTrack.noEndVertex()) // || myTrack.genpartIndex()>=0) + { + HGCal_helpers::Coordinates propcoords; + bool reachesHGCal = toHGCalPropagator.propagate( + myTrack.momentum(), myTrack.vertex().position(), myTrack.charge(), propcoords); + vtx = propcoords.toVector(); + + if (reachesHGCal && vtx.Rho() < eeOuterRadius && vtx.Rho() > eeInnerRadius) { + reachedEE = ReachHGCal::onEESurface; + double dpt = 0; + + for (int i = 0; i < myTrack.nDaughters(); ++i) dpt += myTrack.daughter(i).momentum().pt(); + if (abs(myTrack.type()) == 11) fbrem = dpt / myTrack.momentum().pt(); + } else if (reachesHGCal && vtx.Rho() > eeOuterRadius) + reachedEE = ReachHGCal::outsideEESurface; + + HGCal_helpers::SimpleTrackPropagator indiv_particleProp(aField_); + for (unsigned il = 1; il <= nlayers; ++il) { + const float charge = myTrack.charge(); + indiv_particleProp.setPropagationTargetZ(triggerTools_.getLayerZ(il)); + HGCal_helpers::Coordinates propCoords; + indiv_particleProp.propagate(myTrack.momentum(), myTrack.vertex().position(), charge, + propCoords); + + xp.push_back(propCoords.x); + yp.push_back(propCoords.y); + zp.push_back(propCoords.z); + } + } else { + vtx = myTrack.endVertex().position(); + } + auto orig_vtx = myTrack.vertex().position(); + + allselectedgentracks.push_back(&mySimEvent_->track(i)); + // fill branches + genpart_eta_.push_back(myTrack.momentum().eta()); + genpart_phi_.push_back(myTrack.momentum().phi()); + genpart_pt_.push_back(myTrack.momentum().pt()); + genpart_energy_.push_back(myTrack.momentum().energy()); + genpart_dvx_.push_back(vtx.x()); + genpart_dvy_.push_back(vtx.y()); + genpart_dvz_.push_back(vtx.z()); + + genpart_ovx_.push_back(orig_vtx.x()); + genpart_ovy_.push_back(orig_vtx.y()); + genpart_ovz_.push_back(orig_vtx.z()); + + HGCal_helpers::Coordinates hitsHGCal; + toHGCalPropagator.propagate(myTrack.momentum(), orig_vtx, myTrack.charge(), hitsHGCal); + + genpart_exphi_.push_back(hitsHGCal.phi); + genpart_exeta_.push_back(hitsHGCal.eta); + genpart_exx_.push_back(hitsHGCal.x); + genpart_exy_.push_back(hitsHGCal.y); + + genpart_fbrem_.push_back(fbrem); + genpart_pid_.push_back(myTrack.type()); + genpart_gen_.push_back(myTrack.genpartIndex()); + genpart_reachedEE_.push_back(reachedEE); + genpart_fromBeamPipe_.push_back(true); + + genpart_posx_.push_back(xp); + genpart_posy_.push_back(yp); + genpart_posz_.push_back(zp); + } + + + edm::Handle> genParticlesHandle; + iEvent.getByToken(gen_token_, genParticlesHandle); + gen_n_ = genParticlesHandle->size(); + + for (std::vector::const_iterator it_p = genParticlesHandle->begin(); + it_p != genParticlesHandle->end(); ++it_p) { + gen_eta_.push_back(it_p->eta()); + gen_phi_.push_back(it_p->phi()); + gen_pt_.push_back(it_p->pt()); + gen_energy_.push_back(it_p->energy()); + gen_charge_.push_back(it_p->charge()); + gen_pdgid_.push_back(it_p->pdgId()); + gen_status_.push_back(it_p->status()); + std::vector daughters(it_p->daughterRefVector().size(), 0); + for (unsigned j = 0; j < it_p->daughterRefVector().size(); ++j) { + daughters[j] = static_cast(it_p->daughterRefVector().at(j).key()); + } + gen_daughters_.push_back(daughters); + } + + + // associate gen particles to mothers + genpart_mother_.resize(genpart_posz_.size(), -1); + for (size_t i = 0; i < allselectedgentracks.size(); i++) { + const auto tracki = allselectedgentracks.at(i); + + for (size_t j = i + 1; j < allselectedgentracks.size(); j++) { + const auto trackj = allselectedgentracks.at(j); + + if (!tracki->noMother()) { + if (&tracki->mother() == trackj) genpart_mother_.at(i) = j; + } + if (!trackj->noMother()) { + if (&trackj->mother() == tracki) genpart_mother_.at(j) = i; + } + } + } + + for(const auto& PVI : PupInfo) { - gen_id_.emplace_back(particle.pdgId()); - gen_status_.emplace_back(particle.status()); - gen_energy_.emplace_back(particle.energy()); - gen_pt_.emplace_back(particle.pt()); - gen_eta_.emplace_back(particle.eta()); - gen_phi_.emplace_back(particle.phi()); + if(PVI.getBunchCrossing() == 0) + { + gen_PUNumInt_ = PVI.getPU_NumInteractions(); + gen_TrueNumInt_ = PVI.getTrueNumInteractions(); + } } } @@ -84,15 +454,50 @@ void HGCalTriggerNtupleGen:: clear() { + gen_n_ = 0; - gen_id_.clear(); - gen_status_.clear(); - gen_energy_.clear(); - gen_pt_.clear(); - gen_eta_.clear(); - gen_phi_.clear(); -} + gen_PUNumInt_ = 0; + gen_TrueNumInt_ = 0.; + vtx_x_ = 0; + vtx_y_ = 0; + vtx_z_ = 0; + // + genpart_eta_.clear(); + genpart_phi_.clear(); + genpart_pt_.clear(); + genpart_energy_.clear(); + genpart_dvx_.clear(); + genpart_dvy_.clear(); + genpart_dvz_.clear(); + genpart_ovx_.clear(); + genpart_ovy_.clear(); + genpart_ovz_.clear(); + genpart_exx_.clear(); + genpart_exy_.clear(); + genpart_mother_.clear(); + genpart_exphi_.clear(); + genpart_exeta_.clear(); + genpart_fbrem_.clear(); + genpart_pid_.clear(); + genpart_gen_.clear(); + genpart_reachedEE_.clear(); + genpart_fromBeamPipe_.clear(); + genpart_posx_.clear(); + genpart_posy_.clear(); + genpart_posz_.clear(); + //////////////////// + // reco::GenParticles + // + gen_eta_.clear(); + gen_phi_.clear(); + gen_pt_.clear(); + gen_energy_.clear(); + gen_charge_.clear(); + gen_pdgid_.clear(); + gen_status_.clear(); + gen_daughters_.clear(); +} diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenTau.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenTau.cc index 364db91486dfd..206bfcbdaab50 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenTau.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleGenTau.cc @@ -22,10 +22,13 @@ class HGCalTriggerNtupleGenTau : public HGCalTriggerNtupleBase bool isStableLepton( const reco::GenParticle & daughter ) const; bool isElectron( const reco::GenParticle & daughter ) const; bool isMuon( const reco::GenParticle & daughter ) const; - bool isChargedPion( const reco::GenParticle & daughter ) const; + bool isChargedHadron( const reco::GenParticle & daughter ) const; + bool isChargedHadronFromResonance( const reco::GenParticle & daughter ) const; bool isNeutralPion( const reco::GenParticle & daughter ) const; + bool isNeutralPionFromResonance( const reco::GenParticle & daughter ) const; bool isIntermediateResonance( const reco::GenParticle & daughter ) const; bool isGamma( const reco::GenParticle & daughter ) const; + bool isStableNeutralHadron( const reco::GenParticle & daughter ) const; edm::EDGetToken gen_token_; bool isPythia8generator_; @@ -44,6 +47,7 @@ class HGCalTriggerNtupleGenTau : public HGCalTriggerNtupleBase std::vector gentau_decayMode_; std::vector gentau_totNproducts_; std::vector gentau_totNgamma_; + std::vector gentau_totNpiZero_; std::vector gentau_totNcharged_; std::vector > gentau_products_pt_; @@ -92,6 +96,7 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("gentau_decayMode", &gentau_decayMode_); tree.Branch("gentau_totNproducts", &gentau_totNproducts_); tree.Branch("gentau_totNgamma", &gentau_totNgamma_); + tree.Branch("gentau_totNpiZero", &gentau_totNpiZero_); tree.Branch("gentau_totNcharged", &gentau_totNcharged_); } @@ -101,11 +106,16 @@ bool HGCalTriggerNtupleGenTau::isGoodTau( const reco::GenParticle& candidate ) c } -bool HGCalTriggerNtupleGenTau::isChargedPion( const reco::GenParticle& candidate ) const { - return ( std::abs(candidate.pdgId()) == 211 && candidate.status()==1 +bool HGCalTriggerNtupleGenTau::isChargedHadron( const reco::GenParticle& candidate ) const { + return ( (std::abs(candidate.pdgId()) == 211 || std::abs(candidate.pdgId()) == 321 ) && candidate.status()==1 && candidate.isDirectPromptTauDecayProductFinalState() && candidate.isLastCopy() ); } +bool HGCalTriggerNtupleGenTau::isChargedHadronFromResonance( const reco::GenParticle& candidate ) const { + return ( (std::abs(candidate.pdgId()) == 211 || std::abs(candidate.pdgId()) == 321 ) && candidate.status()==1 + && candidate.isLastCopy() ); +} + bool HGCalTriggerNtupleGenTau::isStableLepton( const reco::GenParticle& candidate ) const { @@ -133,17 +143,29 @@ bool HGCalTriggerNtupleGenTau::isNeutralPion( const reco::GenParticle& candidate } +bool HGCalTriggerNtupleGenTau::isNeutralPionFromResonance( const reco::GenParticle& candidate ) const +{ + return ( std::abs(candidate.pdgId()) == 111 && candidate.status()==2 && candidate.statusFlags().isTauDecayProduct() ); +} + + bool HGCalTriggerNtupleGenTau::isGamma( const reco::GenParticle& candidate ) const { return ( std::abs(candidate.pdgId()) == 22 && candidate.status()==1 && candidate.statusFlags().isTauDecayProduct() && !candidate.isDirectPromptTauDecayProductFinalState() && candidate.isLastCopy() ); } - bool HGCalTriggerNtupleGenTau::isIntermediateResonance( const reco::GenParticle& candidate ) const { return ( ( std::abs(candidate.pdgId()) == 213 || std::abs(candidate.pdgId()) == 20213 || std::abs(candidate.pdgId()) == 24 ) - && candidate.isDirectPromptTauDecayProductFinalState() && candidate.status() == 2 ); + && candidate.status() == 2 ); +} + + +bool HGCalTriggerNtupleGenTau::isStableNeutralHadron( const reco::GenParticle& candidate ) const +{ + return ( !( std::abs(candidate.pdgId())>10 && std::abs(candidate.pdgId())<17) && !isChargedHadron(candidate) + && candidate.isDirectPromptTauDecayProductFinalState() && candidate.status() == 1 ); } @@ -159,10 +181,8 @@ fill(const edm::Event& e, const edm::EventSetup& es) for(const auto& particle : gen_particles) { - /* select good taus */ if( isGoodTau( particle ) ){ - LorentzVector tau_p4vis(0.,0.,0.,0.); gentau_pt_.emplace_back( particle.pt() ); gentau_eta_.emplace_back( particle.eta() ); @@ -201,13 +221,13 @@ fill(const edm::Event& e, const edm::EventSetup& es) finalProds.push_back( daughter ); } - if( isChargedPion( *daughter ) ){ + else if( isChargedHadron( *daughter ) ){ n_pi++; tau_p4vis+=(daughter->p4()); finalProds.push_back( daughter ); } - if( isNeutralPion( *daughter ) ){ + else if( isNeutralPion( *daughter ) ){ n_piZero++; const reco::GenParticleRefVector& grandaughters = daughter->daughterRefVector(); for( const auto& grandaughter : grandaughters ){ @@ -219,17 +239,21 @@ fill(const edm::Event& e, const edm::EventSetup& es) } } + else if( isStableNeutralHadron( *daughter ) ){ + tau_p4vis+=(daughter->p4()); + finalProds.push_back( daughter ); + } /* Here the selection of the decay product according to the Pythia6 decayTree */ if( !isPythia8generator_ ){ if( isIntermediateResonance( *daughter ) ){ const reco::GenParticleRefVector& grandaughters = daughter->daughterRefVector(); for( const auto& grandaughter : grandaughters ){ - if( isChargedPion( *grandaughter ) ){ + if( isChargedHadron( *grandaughter ) || isChargedHadronFromResonance( *grandaughter ) ){ n_pi++; tau_p4vis+=(grandaughter->p4()); finalProds.push_back( daughter ); } - if( isNeutralPion( *grandaughter ) ){ + else if( isNeutralPion( *grandaughter ) || isNeutralPionFromResonance( *grandaughter ) ){ n_piZero++; const reco::GenParticleRefVector& descendants = grandaughter->daughterRefVector(); for( const auto& descendant : descendants ){ @@ -264,6 +288,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) gentau_vis_mass_.emplace_back(tau_p4vis.M()); gentau_totNproducts_.emplace_back(n_pi + n_gamma); gentau_totNgamma_.emplace_back(n_gamma); + gentau_totNpiZero_.emplace_back(n_piZero); gentau_totNcharged_.emplace_back(n_pi); gentau_products_pt_.emplace_back(tau_products_pt); diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc index 5ba066a2453e7..9f5f890526240 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCClusters.cc @@ -1,10 +1,12 @@ #include #include "DataFormats/L1THGCal/interface/HGCalCluster.h" +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" @@ -21,16 +23,25 @@ class HGCalTriggerNtupleHGCClusters : public HGCalTriggerNtupleBase virtual void clear() override final; - edm::EDGetToken clusters_token_; + bool filter_clusters_in_multiclusters_; + edm::EDGetToken clusters_token_, multiclusters_token_; + HGCalTriggerTools triggerTools_; + int cl_n_ ; + std::vector cl_id_; + std::vector cl_mipPt_; std::vector cl_pt_; std::vector cl_energy_; std::vector cl_eta_; std::vector cl_phi_; std::vector cl_layer_; - std::vector cl_ncells_; - std::vector> cl_cells_; + std::vector cl_subdet_; + std::vector cl_cells_n_; + std::vector> cl_cells_id_; + std::vector cl_multicluster_id_; + std::vector cl_multicluster_pt_; + }; @@ -42,6 +53,7 @@ DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleHGCClusters:: HGCalTriggerNtupleHGCClusters(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) { + filter_clusters_in_multiclusters_ = conf.getParameter("FilterClustersInMulticlusters"); } void @@ -49,15 +61,21 @@ HGCalTriggerNtupleHGCClusters:: initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector) { clusters_token_ = collector.consumes(conf.getParameter("Clusters")); + multiclusters_token_ = collector.consumes(conf.getParameter("Multiclusters")); tree.Branch("cl_n", &cl_n_, "cl_n/I"); + tree.Branch("cl_id", &cl_id_); + tree.Branch("cl_mipPt", &cl_mipPt_); tree.Branch("cl_pt", &cl_pt_); tree.Branch("cl_energy", &cl_energy_); tree.Branch("cl_eta", &cl_eta_); - tree.Branch("cl_phi", &cl_phi_); + tree.Branch("cl_phi", &cl_phi_); tree.Branch("cl_layer", &cl_layer_); - tree.Branch("cl_ncells", &cl_ncells_); - tree.Branch("cl_cells", &cl_cells_); + tree.Branch("cl_subdet", &cl_subdet_); + tree.Branch("cl_cells_n", &cl_cells_n_); + tree.Branch("cl_cells_id", &cl_cells_id_); + tree.Branch("cl_multicluster_id", &cl_multicluster_id_); + tree.Branch("cl_multicluster_pt", &cl_multicluster_pt_); } void @@ -69,27 +87,55 @@ fill(const edm::Event& e, const edm::EventSetup& es) edm::Handle clusters_h; e.getByToken(clusters_token_, clusters_h); const l1t::HGCalClusterBxCollection& clusters = *clusters_h; + edm::Handle multiclusters_h; + e.getByToken(multiclusters_token_, multiclusters_h); + const l1t::HGCalMulticlusterBxCollection& multiclusters = *multiclusters_h; // retrieve geometry edm::ESHandle geometry; es.get().get(geometry); + triggerTools_.setEventSetup(es); + + // Associate cells to clusters + std::unordered_map cluster2multicluster; + for(auto mcl_itr=multiclusters.begin(0); mcl_itr!=multiclusters.end(0); mcl_itr++) + { + // loop on 2D clusters inside 3D clusters + for(const auto& cl_ptr : mcl_itr->constituents()) + { + cluster2multicluster.emplace(cl_ptr->detId(), mcl_itr); + } + } + + + clear(); for(auto cl_itr=clusters.begin(0); cl_itr!=clusters.end(0); cl_itr++) { + auto mcl_itr = cluster2multicluster.find(cl_itr->detId()); + uint32_t mcl_id = (mcl_itr!=cluster2multicluster.end() ? mcl_itr->second->detId() : 0); + float mcl_pt = (mcl_itr!=cluster2multicluster.end() ? mcl_itr->second->pt() : 0.); + if(filter_clusters_in_multiclusters_ && mcl_id==0) continue; cl_n_++; + cl_mipPt_.emplace_back(cl_itr->mipPt()); // physical values cl_pt_.emplace_back(cl_itr->pt()); cl_energy_.emplace_back(cl_itr->energy()); cl_eta_.emplace_back(cl_itr->eta()); cl_phi_.emplace_back(cl_itr->phi()); - cl_layer_.emplace_back(cl_itr->layer()); - cl_ncells_.emplace_back(cl_itr->constituents().size()); + + cl_id_.emplace_back(cl_itr->detId()); + cl_layer_.emplace_back(triggerTools_.getLayerWithOffset(cl_itr->detId())); + cl_subdet_.emplace_back(cl_itr->subdetId()); + cl_cells_n_.emplace_back(cl_itr->constituents().size()); // Retrieve indices of trigger cells inside cluster - cl_cells_.emplace_back(cl_itr->constituents().size()); + cl_cells_id_.emplace_back(cl_itr->constituents().size()); std::transform(cl_itr->constituents_begin(), cl_itr->constituents_end(), - cl_cells_.back().begin(), [](const edm::Ptr& tc){return tc.key();} + cl_cells_id_.back().begin(), [](const edm::Ptr& tc){return tc->detId();} ); + cl_multicluster_id_.emplace_back(mcl_id); + cl_multicluster_pt_.emplace_back(mcl_pt); } } @@ -99,15 +145,16 @@ HGCalTriggerNtupleHGCClusters:: clear() { cl_n_ = 0; + cl_id_.clear(); + cl_mipPt_.clear(); cl_pt_.clear(); cl_energy_.clear(); cl_eta_.clear(); cl_phi_.clear(); cl_layer_.clear(); - cl_ncells_.clear(); - cl_cells_.clear(); + cl_subdet_.clear(); + cl_cells_n_.clear(); + cl_cells_id_.clear(); + cl_multicluster_id_.clear(); + cl_multicluster_pt_.clear(); } - - - - diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc index c6e98cb712c85..e71c5ae8a4adf 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCDigis.cc @@ -13,6 +13,8 @@ #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" + class HGCalTriggerNtupleHGCDigis : public HGCalTriggerNtupleBase { @@ -31,6 +33,8 @@ class HGCalTriggerNtupleHGCDigis : public HGCalTriggerNtupleBase bool is_Simhit_comp_; edm::EDGetToken SimHits_inputee_, SimHits_inputfh_, SimHits_inputbh_; + HGCalTriggerTools triggerTools_; + int hgcdigi_n_ ; std::vector hgcdigi_id_; std::vector hgcdigi_subdet_; @@ -60,7 +64,7 @@ class HGCalTriggerNtupleHGCDigis : public HGCalTriggerNtupleBase std::vector bhdigi_simenergy_; edm::ESHandle triggerGeometry_; - + }; DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, @@ -80,7 +84,7 @@ HGCalTriggerNtupleHGCDigis:: initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector) { - ee_token_ = collector.consumes(conf.getParameter("HGCDigisEE")); + ee_token_ = collector.consumes(conf.getParameter("HGCDigisEE")); fh_token_ = collector.consumes(conf.getParameter("HGCDigisFH")); bh_token_ = collector.consumes(conf.getParameter("HGCDigisBH")); if (is_Simhit_comp_) { @@ -95,10 +99,10 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("hgcdigi_layer", &hgcdigi_layer_); tree.Branch("hgcdigi_wafer", &hgcdigi_wafer_); tree.Branch("hgcdigi_wafertype", &hgcdigi_wafertype_); - tree.Branch("hgcdigi_cell", &hgcdigi_cell_); - tree.Branch("hgcdigi_eta", &hgcdigi_eta_); - tree.Branch("hgcdigi_phi", &hgcdigi_phi_); - tree.Branch("hgcdigi_z", &hgcdigi_z_); + tree.Branch("hgcdigi_cell", &hgcdigi_cell_); + tree.Branch("hgcdigi_eta", &hgcdigi_eta_); + tree.Branch("hgcdigi_phi", &hgcdigi_phi_); + tree.Branch("hgcdigi_z", &hgcdigi_z_); tree.Branch("hgcdigi_data", &hgcdigi_data_); tree.Branch("hgcdigi_isadc", &hgcdigi_isadc_); if (is_Simhit_comp_) tree.Branch("hgcdigi_simenergy", &hgcdigi_simenergy_); @@ -109,10 +113,10 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("bhdigi_zside", &bhdigi_side_); tree.Branch("bhdigi_layer", &bhdigi_layer_); tree.Branch("bhdigi_ieta", &bhdigi_ieta_); - tree.Branch("bhdigi_iphi", &bhdigi_iphi_); - tree.Branch("bhdigi_eta", &bhdigi_eta_); - tree.Branch("bhdigi_phi", &bhdigi_phi_); - tree.Branch("bhdigi_z", &bhdigi_z_); + tree.Branch("bhdigi_iphi", &bhdigi_iphi_); + tree.Branch("bhdigi_eta", &bhdigi_eta_); + tree.Branch("bhdigi_phi", &bhdigi_phi_); + tree.Branch("bhdigi_z", &bhdigi_z_); tree.Branch("bhdigi_data", &bhdigi_data_); if (is_Simhit_comp_) tree.Branch("bhdigi_simenergy", &bhdigi_simenergy_); } @@ -133,12 +137,14 @@ fill(const edm::Event& e, const edm::EventSetup& es) e.getByToken(bh_token_, bh_digis_h); const HGCBHDigiCollection& bh_digis = *bh_digis_h; + triggerTools_.setEventSetup(es); + // sim hit association std::unordered_map simhits_ee; - std::unordered_map simhits_fh; - std::unordered_map simhits_bh; + std::unordered_map simhits_fh; + std::unordered_map simhits_bh; if (is_Simhit_comp_) simhits(e, simhits_ee, simhits_fh, simhits_bh); - + clear(); hgcdigi_n_ = ee_digis.size() + fh_digis.size(); hgcdigi_id_.reserve(hgcdigi_n_); @@ -166,7 +172,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) bhdigi_phi_.reserve(bhdigi_n_); bhdigi_z_.reserve(bhdigi_n_); if (is_Simhit_comp_) bhdigi_simenergy_.reserve(bhdigi_n_); - + const int kIntimeSample = 2; for(const auto& digi : ee_digis) { @@ -174,7 +180,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) hgcdigi_id_.emplace_back(id.rawId()); hgcdigi_subdet_.emplace_back(ForwardSubdetector::HGCEE); hgcdigi_side_.emplace_back(id.zside()); - hgcdigi_layer_.emplace_back(id.layer()); + hgcdigi_layer_.emplace_back(triggerTools_.getLayerWithOffset(id)); hgcdigi_wafer_.emplace_back(id.wafer()); hgcdigi_wafertype_.emplace_back(id.waferType()); hgcdigi_cell_.emplace_back(id.cell()); @@ -182,7 +188,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) hgcdigi_eta_.emplace_back(cellpos.eta()); hgcdigi_phi_.emplace_back(cellpos.phi()); hgcdigi_z_.emplace_back(cellpos.z()); - hgcdigi_data_.emplace_back(digi[kIntimeSample].data()); + hgcdigi_data_.emplace_back(digi[kIntimeSample].data()); int is_adc=0; if (!(digi[kIntimeSample].mode())) is_adc =1; hgcdigi_isadc_.emplace_back(is_adc); @@ -190,17 +196,17 @@ fill(const edm::Event& e, const edm::EventSetup& es) double hit_energy=0; auto itr = simhits_ee.find(id); if(itr!=simhits_ee.end())hit_energy = itr->second; - hgcdigi_simenergy_.emplace_back(hit_energy); + hgcdigi_simenergy_.emplace_back(hit_energy); } } - + for(const auto& digi : fh_digis) { const HGCalDetId id(digi.id()); hgcdigi_id_.emplace_back(id.rawId()); hgcdigi_subdet_.emplace_back(ForwardSubdetector::HGCHEF); hgcdigi_side_.emplace_back(id.zside()); - hgcdigi_layer_.emplace_back(id.layer()); + hgcdigi_layer_.emplace_back(triggerTools_.getLayerWithOffset(id)); hgcdigi_wafer_.emplace_back(id.wafer()); hgcdigi_wafertype_.emplace_back(id.waferType()); hgcdigi_cell_.emplace_back(id.cell()); @@ -208,7 +214,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) hgcdigi_eta_.emplace_back(cellpos.eta()); hgcdigi_phi_.emplace_back(cellpos.phi()); hgcdigi_z_.emplace_back(cellpos.z()); - hgcdigi_data_.emplace_back(digi[kIntimeSample].data()); + hgcdigi_data_.emplace_back(digi[kIntimeSample].data()); int is_adc=0; if (!(digi[kIntimeSample].mode())) is_adc =1; hgcdigi_isadc_.emplace_back(is_adc); @@ -216,7 +222,7 @@ fill(const edm::Event& e, const edm::EventSetup& es) double hit_energy=0; auto itr = simhits_fh.find(id); if(itr!=simhits_fh.end())hit_energy = itr->second; - hgcdigi_simenergy_.emplace_back(hit_energy); + hgcdigi_simenergy_.emplace_back(hit_energy); } } @@ -226,19 +232,19 @@ fill(const edm::Event& e, const edm::EventSetup& es) bhdigi_id_.emplace_back(id.rawId()); bhdigi_subdet_.emplace_back(id.subdetId()); bhdigi_side_.emplace_back(id.zside()); - bhdigi_layer_.emplace_back(id.depth()); + bhdigi_layer_.emplace_back(triggerTools_.getLayerWithOffset(id)); bhdigi_ieta_.emplace_back(id.ieta()); bhdigi_iphi_.emplace_back(id.iphi()); GlobalPoint cellpos = triggerGeometry_->bhGeometry().getPosition(id.rawId()); bhdigi_eta_.emplace_back(cellpos.eta()); bhdigi_phi_.emplace_back(cellpos.phi()); bhdigi_z_.emplace_back(cellpos.z()); - bhdigi_data_.emplace_back(digi[kIntimeSample].data()); + bhdigi_data_.emplace_back(digi[kIntimeSample].data()); if (is_Simhit_comp_) { double hit_energy=0; auto itr = simhits_bh.find(id); if(itr!=simhits_bh.end())hit_energy = itr->second; - bhdigi_simenergy_.emplace_back(hit_energy); + bhdigi_simenergy_.emplace_back(hit_energy); } } } @@ -257,13 +263,13 @@ simhits(const edm::Event& e, std::unordered_map& simhits_ee, s edm::Handle bh_simhits_h; e.getByToken(SimHits_inputbh_,bh_simhits_h); const edm::PCaloHitContainer& bh_simhits = *bh_simhits_h; - + //EE int layer=0,cell=0, sec=0, subsec=0, zp=0,subdet=0; ForwardSubdetector mysubdet; - - for( const auto& simhit : ee_simhits ) { - HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); + + for( const auto& simhit : ee_simhits ) { + HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); mysubdet = (ForwardSubdetector)(subdet); std::pair recoLayerCell = triggerGeometry_->eeTopology().dddConstants().simToReco(cell,layer,sec,triggerGeometry_->eeTopology().detectorType()); cell = recoLayerCell.first; @@ -277,9 +283,9 @@ simhits(const edm::Event& e, std::unordered_map& simhits_ee, s // FH layer=0; cell=0; sec=0; subsec=0; zp=0; subdet=0; - - for( const auto& simhit : fh_simhits ) { - HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); + + for( const auto& simhit : fh_simhits ) { + HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); mysubdet = (ForwardSubdetector)(subdet); std::pair recoLayerCell = triggerGeometry_->fhTopology().dddConstants().simToReco(cell,layer,sec,triggerGeometry_->fhTopology().detectorType()); cell = recoLayerCell.first; @@ -289,14 +295,14 @@ simhits(const edm::Event& e, std::unordered_map& simhits_ee, s } auto itr_insert = simhits_fh.emplace(HGCalDetId(mysubdet,zp,layer,subsec,sec,cell), 0.); itr_insert.first->second += simhit.energy(); - } + } // BH - for( const auto& simhit : bh_simhits ) { + for( const auto& simhit : bh_simhits ) { HcalDetId id = HcalHitRelabeller::relabel(simhit.id(), triggerGeometry_->bhTopology().dddConstants()); if (id.subdetId()!=HcalEndcap) continue; auto itr_insert = simhits_bh.emplace(id, 0.); itr_insert.first->second += simhit.energy(); - } + } } @@ -332,7 +338,3 @@ clear() bhdigi_data_.clear(); if (is_Simhit_comp_) bhdigi_simenergy_.clear(); } - - - - diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc index 775c070810e74..524da02bfa3bb 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc @@ -1,4 +1,3 @@ - #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" @@ -19,16 +18,30 @@ class HGCalTriggerNtupleHGCMulticlusters : public HGCalTriggerNtupleBase private: virtual void clear() override final; - edm::EDGetToken multiclusters_token_; int cl3d_n_ ; + std::vector cl3d_id_; std::vector cl3d_pt_; std::vector cl3d_energy_; std::vector cl3d_eta_; std::vector cl3d_phi_; - std::vector cl3d_nclu_; - std::vector> cl3d_clusters_; + std::vector cl3d_clusters_n_; + std::vector> cl3d_clusters_id_; + // cluster shower shapes + std::vector cl3d_showerlength_; + std::vector cl3d_coreshowerlength_; + std::vector cl3d_firstlayer_; + std::vector cl3d_maxlayer_; + std::vector cl3d_seetot_; + std::vector cl3d_seemax_; + std::vector cl3d_spptot_; + std::vector cl3d_sppmax_; + std::vector cl3d_szz_; + std::vector cl3d_srrtot_; + std::vector cl3d_srrmax_; + std::vector cl3d_srrmean_; + std::vector cl3d_emaxe_; }; DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, @@ -48,12 +61,26 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& multiclusters_token_ = collector.consumes(conf.getParameter("Multiclusters")); tree.Branch("cl3d_n", &cl3d_n_, "cl3d_n/I"); + tree.Branch("cl3d_id", &cl3d_id_); tree.Branch("cl3d_pt", &cl3d_pt_); tree.Branch("cl3d_energy", &cl3d_energy_); tree.Branch("cl3d_eta", &cl3d_eta_); tree.Branch("cl3d_phi", &cl3d_phi_); - tree.Branch("cl3d_nclu", &cl3d_nclu_); - tree.Branch("cl3d_clusters", &cl3d_clusters_); + tree.Branch("cl3d_clusters_n", &cl3d_clusters_n_); + tree.Branch("cl3d_clusters_id", &cl3d_clusters_id_); + tree.Branch("cl3d_showerlength", &cl3d_showerlength_); + tree.Branch("cl3d_coreshowerlength", &cl3d_coreshowerlength_); + tree.Branch("cl3d_firstlayer", &cl3d_firstlayer_); + tree.Branch("cl3d_maxlayer", &cl3d_maxlayer_); + tree.Branch("cl3d_seetot", &cl3d_seetot_); + tree.Branch("cl3d_seemax", &cl3d_seemax_); + tree.Branch("cl3d_spptot", &cl3d_spptot_); + tree.Branch("cl3d_sppmax", &cl3d_sppmax_); + tree.Branch("cl3d_szz", &cl3d_szz_); + tree.Branch("cl3d_srrtot", &cl3d_srrtot_); + tree.Branch("cl3d_srrmax", &cl3d_srrmax_); + tree.Branch("cl3d_srrmean", &cl3d_srrmean_); + tree.Branch("cl3d_emaxe", &cl3d_emaxe_); } @@ -75,16 +102,31 @@ fill(const edm::Event& e, const edm::EventSetup& es) for(auto cl3d_itr=multiclusters.begin(0); cl3d_itr!=multiclusters.end(0); cl3d_itr++) { cl3d_n_++; + cl3d_id_.emplace_back(cl3d_itr->detId()); // physical values cl3d_pt_.emplace_back(cl3d_itr->pt()); cl3d_energy_.emplace_back(cl3d_itr->energy()); cl3d_eta_.emplace_back(cl3d_itr->eta()); cl3d_phi_.emplace_back(cl3d_itr->phi()); - cl3d_nclu_.emplace_back(cl3d_itr->constituents().size()); + cl3d_clusters_n_.emplace_back(cl3d_itr->constituents().size()); + cl3d_showerlength_.emplace_back(cl3d_itr->showerLength()); + cl3d_coreshowerlength_.emplace_back(cl3d_itr->coreShowerLength()); + cl3d_firstlayer_.emplace_back(cl3d_itr->firstLayer()); + cl3d_maxlayer_.emplace_back(cl3d_itr->maxLayer()); + cl3d_seetot_.emplace_back(cl3d_itr->sigmaEtaEtaTot()); + cl3d_seemax_.emplace_back(cl3d_itr->sigmaEtaEtaMax()); + cl3d_spptot_.emplace_back(cl3d_itr->sigmaPhiPhiTot()); + cl3d_sppmax_.emplace_back(cl3d_itr->sigmaPhiPhiMax()); + cl3d_szz_.emplace_back(cl3d_itr->sigmaZZ()); + cl3d_srrtot_.emplace_back(cl3d_itr->sigmaRRTot()); + cl3d_srrmax_.emplace_back(cl3d_itr->sigmaRRMax()); + cl3d_srrmean_.emplace_back(cl3d_itr->sigmaRRMean()); + cl3d_emaxe_.emplace_back(cl3d_itr->eMax()/cl3d_itr->energy()); + // Retrieve indices of trigger cells inside cluster - cl3d_clusters_.emplace_back(cl3d_itr->constituents().size()); + cl3d_clusters_id_.emplace_back(cl3d_itr->constituents().size()); std::transform(cl3d_itr->constituents_begin(), cl3d_itr->constituents_end(), - cl3d_clusters_.back().begin(), [](const edm::Ptr& cl){return cl.key();} + cl3d_clusters_id_.back().begin(), [](const edm::Ptr& cl){return cl->detId();} ); } } @@ -95,13 +137,26 @@ HGCalTriggerNtupleHGCMulticlusters:: clear() { cl3d_n_ = 0; + cl3d_id_.clear(); cl3d_pt_.clear(); cl3d_energy_.clear(); cl3d_eta_.clear(); cl3d_phi_.clear(); - cl3d_nclu_.clear(); - cl3d_clusters_.clear(); - + cl3d_clusters_n_.clear(); + cl3d_clusters_id_.clear(); + cl3d_showerlength_.clear(); + cl3d_coreshowerlength_.clear(); + cl3d_firstlayer_.clear(); + cl3d_maxlayer_.clear(); + cl3d_seetot_.clear(); + cl3d_seemax_.clear(); + cl3d_spptot_.clear(); + cl3d_sppmax_.clear(); + cl3d_szz_.clear(); + cl3d_srrtot_.clear(); + cl3d_srrmax_.clear(); + cl3d_srrmean_.clear(); + cl3d_emaxe_.clear(); } diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc new file mode 100644 index 0000000000000..d0ac617ada1a1 --- /dev/null +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCPanels.cc @@ -0,0 +1,171 @@ + +#include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h" +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" + + + +class HGCalTriggerNtupleHGCPanels : public HGCalTriggerNtupleBase +{ + + public: + HGCalTriggerNtupleHGCPanels(const edm::ParameterSet& conf); + ~HGCalTriggerNtupleHGCPanels(){}; + virtual void initialize(TTree&, const edm::ParameterSet&, edm::ConsumesCollector&&) override final; + virtual void fill(const edm::Event& e, const edm::EventSetup& es) override final; + + private: + virtual void clear() override final; + + + edm::EDGetToken trigger_cells_token_; + edm::ESHandle geometry_; + + int panel_n_ ; + std::vector panel_id_; + std::vector panel_zside_; + std::vector panel_layer_; + std::vector panel_sector_; + std::vector panel_number_; + std::vector panel_tc_n_; + std::vector > panel_tc_id_; + std::vector > panel_tc_mod_; + std::vector > panel_tc_third_; + std::vector > panel_tc_cell_; + std::vector > panel_tc_mipPt_; + std::vector > panel_tc_pt_; + + private: + static const unsigned kPanel_offset_ = 0; + static const unsigned kPanel_mask_ = 0x1F; + static const unsigned kSector_offset_ = 5; + static const unsigned kSector_mask_ = 0x7; + static const unsigned kThird_offset_ = 4; + static const unsigned kThird_mask_ = 0x3; + static const unsigned kCell_mask_ = 0xF; + +}; + +DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, + HGCalTriggerNtupleHGCPanels, + "HGCalTriggerNtupleHGCPanels" ); + + +HGCalTriggerNtupleHGCPanels:: +HGCalTriggerNtupleHGCPanels(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) +{ +} + +void +HGCalTriggerNtupleHGCPanels:: +initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector) +{ + trigger_cells_token_ = collector.consumes(conf.getParameter("TriggerCells")); + + tree.Branch("panel_n", &panel_n_, "panel_n/I"); + tree.Branch("panel_id", &panel_id_); + tree.Branch("panel_zside", &panel_zside_); + tree.Branch("panel_layer", &panel_layer_); + tree.Branch("panel_sector", &panel_sector_); + tree.Branch("panel_number", &panel_number_); + tree.Branch("panel_tc_n", &panel_tc_n_); + tree.Branch("panel_tc_id", &panel_tc_id_); + tree.Branch("panel_tc_mod", &panel_tc_mod_); + tree.Branch("panel_tc_third", &panel_tc_third_); + tree.Branch("panel_tc_cell", &panel_tc_cell_); + tree.Branch("panel_tc_mipPt", &panel_tc_mipPt_); + tree.Branch("panel_tc_pt", &panel_tc_pt_); + +} + +void +HGCalTriggerNtupleHGCPanels:: +fill(const edm::Event& e, const edm::EventSetup& es) +{ + + // retrieve trigger cells + edm::Handle trigger_cells_h; + e.getByToken(trigger_cells_token_, trigger_cells_h); + const l1t::HGCalTriggerCellBxCollection& trigger_cells = *trigger_cells_h; + + // retrieve geometry + es.get().get(geometry_); + + clear(); + + // Regroup trigger cells by panel + std::unordered_map< uint32_t,vector > panelids_tcs; + for(auto tc_itr=trigger_cells.begin(0); tc_itr!=trigger_cells.end(0); tc_itr++) + { + if(tc_itr->hwPt()>0 && tc_itr->subdetId()!=ForwardSubdetector::HGCHEB) + { + HGCalDetId id(tc_itr->detId()); + HGCalDetId panelid(geometry_->getModuleFromTriggerCell(id)); + panelids_tcs[panelid].push_back(tc_itr); + } + } + for (const auto& panelid_tcs : panelids_tcs) + { + panel_n_++; + HGCalDetId panelid(panelid_tcs.first); + int panel_sector = (panelid.wafer()>>kSector_offset_) & kSector_mask_ ; + int panel_number = (panelid.wafer()>>kPanel_offset_) & kPanel_mask_ ; + const auto& tcs = panelid_tcs.second; + panel_id_.emplace_back(panelid); + panel_zside_.emplace_back(panelid.zside()); + panel_layer_.emplace_back(geometry_->triggerLayer(panelid)); + panel_sector_.emplace_back(panel_sector); + panel_number_.emplace_back(panel_number); + panel_tc_n_.emplace_back(tcs.size()); + panel_tc_id_.emplace_back(); + panel_tc_mod_.emplace_back(); + panel_tc_third_.emplace_back(); + panel_tc_cell_.emplace_back(); + panel_tc_mipPt_.emplace_back(); + panel_tc_pt_.emplace_back(); + + std::unordered_map modules_mipPt; + std::unordered_map thirds_mipPt; + for (const auto& tc : tcs) + { + panel_tc_id_.back().push_back(tc->detId()); + panel_tc_mipPt_.back().push_back(tc->mipPt()); + panel_tc_pt_.back().push_back(tc->pt()); + HGCalDetId tc_detid(tc->detId()); + unsigned module_id = tc_detid.wafer(); + unsigned third_id = (tc_detid.cell()>>kThird_offset_) & kThird_mask_; + unsigned cell_id = tc_detid.cell() & kCell_mask_; + panel_tc_mod_.back().push_back(module_id); + panel_tc_third_.back().push_back(third_id); + panel_tc_cell_.back().push_back(cell_id); + modules_mipPt[module_id] += tc->mipPt(); + thirds_mipPt[third_id] += tc->mipPt(); + } + } +} + + +void +HGCalTriggerNtupleHGCPanels:: +clear() +{ + panel_n_ = 0; + panel_id_.clear(); + panel_zside_.clear(); + panel_layer_.clear(); + panel_sector_.clear(); + panel_number_.clear(); + panel_tc_n_.clear(); + panel_tc_id_.clear(); + panel_tc_mod_.clear(); + panel_tc_third_.clear(); + panel_tc_cell_.clear(); + panel_tc_mipPt_.clear(); + panel_tc_pt_.clear(); +} + + + + diff --git a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc index 73dc9b50e6b4d..36766c0d0e64f 100644 --- a/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc +++ b/L1Trigger/L1THGCal/plugins/ntuples/HGCalTriggerNtupleHGCTriggerCells.cc @@ -1,9 +1,15 @@ +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" +#include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" +#include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h" #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h" +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" #include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" #include "L1Trigger/L1THGCal/interface/HGCalTriggerNtupleBase.h" +#include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" @@ -17,10 +23,22 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase virtual void fill(const edm::Event& e, const edm::EventSetup& es) override final; private: + double calibrate(double, int, int); + void simhits(const edm::Event& e, std::unordered_map& simhits_ee, std::unordered_map& simhits_fh, std::unordered_map& simhits_bh); virtual void clear() override final; + HGCalTriggerTools triggerTools_; + + edm::EDGetToken trigger_cells_token_, multiclusters_token_; + edm::EDGetToken simhits_ee_token_, simhits_fh_token_, simhits_bh_token_; + bool fill_simenergy_; + bool filter_cells_in_multiclusters_; + double keV2fC_; + std::vector fcPerMip_; + std::vector layerWeights_; + std::vector thicknessCorrections_; + edm::ESHandle geometry_; - edm::EDGetToken trigger_cells_token_; int tc_n_ ; std::vector tc_id_; @@ -31,10 +49,18 @@ class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase std::vector tc_wafertype_ ; std::vector tc_cell_; std::vector tc_data_; + std::vector tc_mipPt_; + std::vector tc_pt_; std::vector tc_energy_; + std::vector tc_simenergy_; std::vector tc_eta_; std::vector tc_phi_; + std::vector tc_x_; + std::vector tc_y_; std::vector tc_z_; + std::vector tc_cluster_id_; + std::vector tc_multicluster_id_; + std::vector tc_multicluster_pt_; }; @@ -46,6 +72,12 @@ DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleHGCTriggerCells:: HGCalTriggerNtupleHGCTriggerCells(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf) { + fill_simenergy_ = conf.getParameter("FillSimEnergy"); + filter_cells_in_multiclusters_ = conf.getParameter("FilterCellsInMulticlusters"); + keV2fC_ = conf.getParameter("keV2fC"); + fcPerMip_ = conf.getParameter>("fcPerMip"); + layerWeights_ = conf.getParameter>("layerWeights"); + thicknessCorrections_ = conf.getParameter>("thicknessCorrections"); } void @@ -53,6 +85,14 @@ HGCalTriggerNtupleHGCTriggerCells:: initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector) { trigger_cells_token_ = collector.consumes(conf.getParameter("TriggerCells")); + multiclusters_token_ = collector.consumes(conf.getParameter("Multiclusters")); + + if (fill_simenergy_) + { + simhits_ee_token_ = collector.consumes(conf.getParameter("eeSimHits")); + simhits_fh_token_ = collector.consumes(conf.getParameter("fhSimHits")); + simhits_bh_token_ = collector.consumes(conf.getParameter("bhSimHits")); + } tree.Branch("tc_n", &tc_n_, "tc_n/I"); tree.Branch("tc_id", &tc_id_); @@ -61,12 +101,20 @@ initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& tree.Branch("tc_layer", &tc_layer_); tree.Branch("tc_wafer", &tc_wafer_); tree.Branch("tc_wafertype", &tc_wafertype_); - tree.Branch("tc_cell", &tc_cell_); + tree.Branch("tc_cell", &tc_cell_); tree.Branch("tc_data", &tc_data_); + tree.Branch("tc_pt", &tc_pt_); + tree.Branch("tc_mipPt", &tc_mipPt_); tree.Branch("tc_energy", &tc_energy_); + if(fill_simenergy_) tree.Branch("tc_simenergy", &tc_simenergy_); tree.Branch("tc_eta", &tc_eta_); tree.Branch("tc_phi", &tc_phi_); + tree.Branch("tc_x", &tc_x_); + tree.Branch("tc_y", &tc_y_); tree.Branch("tc_z", &tc_z_); + tree.Branch("tc_cluster_id", &tc_cluster_id_); + tree.Branch("tc_multicluster_id", &tc_multicluster_id_); + tree.Branch("tc_multicluster_pt", &tc_multicluster_pt_); } @@ -80,35 +128,189 @@ fill(const edm::Event& e, const edm::EventSetup& es) e.getByToken(trigger_cells_token_, trigger_cells_h); const l1t::HGCalTriggerCellBxCollection& trigger_cells = *trigger_cells_h; + // retrieve clusters + edm::Handle multiclusters_h; + e.getByToken(multiclusters_token_, multiclusters_h); + const l1t::HGCalMulticlusterBxCollection& multiclusters = *multiclusters_h; + // retrieve geometry - edm::ESHandle geometry; - es.get().get(geometry); + es.get().get(geometry_); + + // sim hit association + std::unordered_map simhits_ee; + std::unordered_map simhits_fh; + std::unordered_map simhits_bh; + if(fill_simenergy_) simhits(e, simhits_ee, simhits_fh, simhits_bh); + + // Associate cells to clusters + std::unordered_map cell2cluster; + std::unordered_map cell2multicluster; + for(auto mcl_itr=multiclusters.begin(0); mcl_itr!=multiclusters.end(0); mcl_itr++) + { + // loop on 2D clusters inside 3D clusters + for(const auto& cl_ptr : mcl_itr->constituents()) + { + // loop on TC inside 2D clusters + for(const auto& tc_ptr : cl_ptr->constituents()) + { + cell2cluster.emplace(tc_ptr->detId(), cl_ptr->detId()); + cell2multicluster.emplace(tc_ptr->detId(), mcl_itr); + } + } + } + + triggerTools_.setEventSetup(es); clear(); for(auto tc_itr=trigger_cells.begin(0); tc_itr!=trigger_cells.end(0); tc_itr++) { if(tc_itr->hwPt()>0) { + auto cl_itr = cell2cluster.find(tc_itr->detId()); + auto mcl_itr = cell2multicluster.find(tc_itr->detId()); + uint32_t cl_id = (cl_itr!=cell2cluster.end() ? cl_itr->second : 0); + uint32_t mcl_id = (mcl_itr!=cell2multicluster.end() ? mcl_itr->second->detId() : 0); + float mcl_pt = (mcl_itr!=cell2multicluster.end() ? mcl_itr->second->pt() : 0.); + // Filter cells not included in a multicluster, if requested + if(filter_cells_in_multiclusters_ && mcl_id==0) continue; tc_n_++; // hardware data HGCalDetId id(tc_itr->detId()); tc_id_.emplace_back(tc_itr->detId()); tc_subdet_.emplace_back(id.subdetId()); tc_side_.emplace_back(id.zside()); - tc_layer_.emplace_back(id.layer()); + tc_layer_.emplace_back(triggerTools_.getLayerWithOffset(id)); tc_wafer_.emplace_back(id.wafer()); tc_wafertype_.emplace_back(id.waferType()); tc_cell_.emplace_back(id.cell()); tc_data_.emplace_back(tc_itr->hwPt()); + tc_mipPt_.emplace_back(tc_itr->mipPt()); // physical values + tc_pt_.emplace_back(tc_itr->pt()); tc_energy_.emplace_back(tc_itr->energy()); tc_eta_.emplace_back(tc_itr->eta()); tc_phi_.emplace_back(tc_itr->phi()); + tc_x_.emplace_back(tc_itr->position().x()); + tc_y_.emplace_back(tc_itr->position().y()); tc_z_.emplace_back(tc_itr->position().z()); + // Links between TC and clusters + tc_cluster_id_.emplace_back(cl_id); + tc_multicluster_id_.emplace_back(mcl_id); + tc_multicluster_pt_.emplace_back(mcl_pt); + + if(fill_simenergy_) + { + double energy = 0; + int subdet = id.subdetId(); + // search for simhit for all the cells inside the trigger cell + for(uint32_t c_id : geometry_->getCellsFromTriggerCell(id)) + { + switch(subdet) + { + case ForwardSubdetector::HGCEE: + { + auto itr = simhits_ee.find(c_id); + if(itr!=simhits_ee.end()) + { + HGCalDetId detid(c_id); + int thickness = geometry_->eeTopology().dddConstants().waferTypeL(detid.wafer())-1; + int layer = detid.layer(); + energy += calibrate(itr->second, thickness, layer); + } + break; + } + case ForwardSubdetector::HGCHEF: + { + auto itr = simhits_fh.find(c_id); + if(itr!=simhits_fh.end()) + { + HGCalDetId detid(c_id); + int thickness = geometry_->fhTopology().dddConstants().waferTypeL(detid.wafer())-1; + int layer = detid.layer(); + energy += calibrate(itr->second, thickness, layer); + } + break; + } + case ForwardSubdetector::HGCHEB: + { + auto itr = simhits_bh.find(c_id); + if(itr!=simhits_bh.end()) energy += itr->second; + break; + } + default: + break; + } + } + tc_simenergy_.emplace_back(energy); + } } } } +double +HGCalTriggerNtupleHGCTriggerCells:: +calibrate(double energy, int thickness, int layer) +{ + double fcPerMip = fcPerMip_[thickness]; + double thicknessCorrection = thicknessCorrections_[thickness]; + double layerWeight = layerWeights_[layer]; + double TeV2GeV = 1.e3; + return energy*keV2fC_/fcPerMip*layerWeight*TeV2GeV/thicknessCorrection; +} + +void +HGCalTriggerNtupleHGCTriggerCells:: +simhits(const edm::Event& e, std::unordered_map& simhits_ee, std::unordered_map& simhits_fh, std::unordered_map& simhits_bh) +{ + edm::Handle ee_simhits_h; + e.getByToken(simhits_ee_token_,ee_simhits_h); + const edm::PCaloHitContainer& ee_simhits = *ee_simhits_h; + edm::Handle fh_simhits_h; + e.getByToken(simhits_fh_token_,fh_simhits_h); + const edm::PCaloHitContainer& fh_simhits = *fh_simhits_h; + edm::Handle bh_simhits_h; + e.getByToken(simhits_bh_token_,bh_simhits_h); + const edm::PCaloHitContainer& bh_simhits = *bh_simhits_h; + + //EE + int layer=0,cell=0, sec=0, subsec=0, zp=0,subdet=0; + ForwardSubdetector mysubdet; + for( const auto& simhit : ee_simhits ) + { + HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); + mysubdet = (ForwardSubdetector)subdet; + std::pair recoLayerCell = geometry_->eeTopology().dddConstants().simToReco(cell,layer,sec,geometry_->eeTopology().detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) continue; + auto itr_insert = simhits_ee.emplace(HGCalDetId(mysubdet,zp,layer,subsec,sec,cell), 0.); + itr_insert.first->second += simhit.energy(); + } + + // FH + layer=0; cell=0; sec=0; subsec=0; zp=0; subdet=0; + for( const auto& simhit : fh_simhits ) + { + HGCalTestNumbering::unpackHexagonIndex(simhit.id(), subdet, zp, layer, sec, subsec, cell); + mysubdet = (ForwardSubdetector)(subdet); + std::pair recoLayerCell = geometry_->fhTopology().dddConstants().simToReco(cell,layer,sec,geometry_->fhTopology().detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) continue; + auto itr_insert = simhits_fh.emplace(HGCalDetId(mysubdet,zp,layer,subsec,sec,cell), 0.); + itr_insert.first->second += simhit.energy(); + } + + // BH + for( const auto& simhit : bh_simhits ) + { + HcalDetId id = HcalHitRelabeller::relabel(simhit.id(), geometry_->bhTopology().dddConstants()); + if (id.subdetId()!=HcalEndcap) continue; + auto itr_insert = simhits_bh.emplace(id, 0.); + itr_insert.first->second += simhit.energy(); + } +} + void HGCalTriggerNtupleHGCTriggerCells:: @@ -123,12 +325,16 @@ clear() tc_wafertype_.clear(); tc_cell_.clear(); tc_data_.clear(); + tc_mipPt_.clear(); + tc_pt_.clear(); tc_energy_.clear(); + tc_simenergy_.clear(); tc_eta_.clear(); tc_phi_.clear(); + tc_x_.clear(); + tc_y_.clear(); tc_z_.clear(); + tc_cluster_id_.clear(); + tc_multicluster_id_.clear(); + tc_multicluster_pt_.clear(); } - - - - diff --git a/L1Trigger/L1THGCal/python/customCalibration.py b/L1Trigger/L1THGCal/python/customCalibration.py new file mode 100644 index 0000000000000..5a80d330e5851 --- /dev/null +++ b/L1Trigger/L1THGCal/python/customCalibration.py @@ -0,0 +1,20 @@ +import FWCore.ParameterSet.Config as cms +import hgcalLayersCalibrationCoefficients_cfi as layercalibparam + + +def custom_cluster_calibration_global(process, + factor=1.084 + ): + parameters_c2d = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].C2d_parameters + parameters_c2d.calibSF_cluster = cms.double(factor) + parameters_c2d.applyLayerCalibration = cms.bool(False) + return process + + +def custom_cluster_calibration_layers(process, + weights=layercalibparam.TrgLayer_weights + ): + parameters_c2d = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].C2d_parameters + parameters_c2d.layerWeights = weights + parameters_c2d.applyLayerCalibration = cms.bool(True) + return process diff --git a/L1Trigger/L1THGCal/python/customClustering.py b/L1Trigger/L1THGCal/python/customClustering.py new file mode 100644 index 0000000000000..1536cbc285581 --- /dev/null +++ b/L1Trigger/L1THGCal/python/customClustering.py @@ -0,0 +1,61 @@ +import FWCore.ParameterSet.Config as cms + + +def custom_2dclustering_distance(process, + distance=6.,# cm + seed_threshold=5.,# MipT + cluster_threshold=2.# MipT + ): + parameters_c2d = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].C2d_parameters + parameters_c2d.seeding_threshold_silicon = cms.double(seed_threshold) + parameters_c2d.seeding_threshold_scintillator = cms.double(seed_threshold) + parameters_c2d.clustering_threshold_silicon = cms.double(cluster_threshold) + parameters_c2d.clustering_threshold_scintillator = cms.double(cluster_threshold) + parameters_c2d.dR_cluster = cms.double(distance) + parameters_c2d.clusterType = cms.string('dRC2d') + return process + +def custom_2dclustering_topological(process, + seed_threshold=5.,# MipT + cluster_threshold=2.# MipT + ): + parameters_c2d = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].C2d_parameters + parameters_c2d.seeding_threshold_silicon = cms.double(seed_threshold) # MipT + parameters_c2d.seeding_threshold_scintillator = cms.double(seed_threshold) # MipT + parameters_c2d.clustering_threshold_silicon = cms.double(cluster_threshold) # MipT + parameters_c2d.clustering_threshold_scintillator = cms.double(cluster_threshold) # MipT + parameters_c2d.clusterType = cms.string('NNC2d') + return process + +def custom_2dclustering_constrainedtopological(process, + distance=6.,# cm + seed_threshold=5.,# MipT + cluster_threshold=2.# MipT + ): + parameters_c2d = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].C2d_parameters + parameters_c2d.seeding_threshold_silicon = cms.double(seed_threshold) # MipT + parameters_c2d.seeding_threshold_scintillator = cms.double(seed_threshold) # MipT + parameters_c2d.clustering_threshold_silicon = cms.double(cluster_threshold) # MipT + parameters_c2d.clustering_threshold_scintillator = cms.double(cluster_threshold) # MipT + parameters_c2d.dR_cluster = cms.double(distance) # cm + parameters_c2d.clusterType = cms.string('dRNNC2d') + return process + +def custom_3dclustering_distance(process, + distance=0.01 + ): + parameters_c3d = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].C3d_parameters + parameters_c3d.dR_multicluster = cms.double(distance) + parameters_c3d.type_multicluster = cms.string('dRC3d') + return process + + +def custom_3dclustering_dbscan(process, + distance=0.005, + min_points=3 + ): + parameters_c3d = process.hgcalTriggerPrimitiveDigiProducer.BEConfiguration.algorithms[0].C3d_parameters + parameters_c3d.dist_dbscan_multicluster = cms.double(distance) + parameters_c3d.minN_dbscan_multicluster = cms.uint32(min_points) + parameters_c3d.type_multicluster = cms.string('DBSCAN3d') + return process diff --git a/L1Trigger/L1THGCal/python/customTriggerGeometry.py b/L1Trigger/L1THGCal/python/customTriggerGeometry.py index 941c2655c70d8..c55267439e83c 100644 --- a/L1Trigger/L1THGCal/python/customTriggerGeometry.py +++ b/L1Trigger/L1THGCal/python/customTriggerGeometry.py @@ -3,7 +3,7 @@ def custom_geometry_ZoltanSplit_V8(process): process.hgcalTriggerGeometryESProducer.TriggerGeometry.TriggerGeometryName = cms.string('HGCalTriggerGeometryHexLayerBasedImp1') process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_8inch_aligned_192_432_V8_0.txt") - process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TModulesMapping = cms.FileInPath("L1Trigger/L1THGCal/data/panel_mapping_60deg_6mod_0.txt") + process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TModulesMapping = cms.FileInPath("L1Trigger/L1THGCal/data/panel_mapping_tdr_0.txt") process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellNeighborsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_8inch_aligned_192_432_0.txt") process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_BH_3x3_30deg_0.txt") process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellNeighborsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_BH_3x3_30deg_0.txt") @@ -16,7 +16,7 @@ def custom_geometry_ZoltanSplit_V8(process): def custom_geometry_ZoltanSplit_V7(process): process.hgcalTriggerGeometryESProducer.TriggerGeometry.TriggerGeometryName = cms.string('HGCalTriggerGeometryHexLayerBasedImp1') process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_8inch_aligned_192_432_V7_0.txt") - process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TModulesMapping = cms.FileInPath("L1Trigger/L1THGCal/data/panel_mapping_60deg_6mod_0.txt") + process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TModulesMapping = cms.FileInPath("L1Trigger/L1THGCal/data/panel_mapping_tdr_0.txt") process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellNeighborsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_8inch_aligned_192_432_0.txt") process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_BH_3x3_30deg_0.txt") process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellNeighborsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_BH_3x3_30deg_0.txt") diff --git a/L1Trigger/L1THGCal/python/hgcalLayersCalibrationCoefficients_cfi.py b/L1Trigger/L1THGCal/python/hgcalLayersCalibrationCoefficients_cfi.py new file mode 100644 index 0000000000000..ebd349289ec37 --- /dev/null +++ b/L1Trigger/L1THGCal/python/hgcalLayersCalibrationCoefficients_cfi.py @@ -0,0 +1,165 @@ +import FWCore.ParameterSet.Config as cms + +AllLayer_weights = cms.vdouble(0.0, + 0.0158115, + 0.0286877, + 0.017707, + 0.00884719, + 0.00921472, + 0.00654193, + 0.00737344, + 0.00737619, + 0.0090818, + 0.00776757, + 0.011098, + 0.00561986, + 0.012015, + 0.0105088, + 0.00834435, + 0.0113901, + 0.00995654, + 0.0120987, + 0.00708785, + 0.0101533, + 0.0108289, + 0.0143815, + 0.0145606, + 0.0133204, + 0.0137476, + 0.00911436, + 0.0275028, + 0.0338628, + 0.0674028, + 0.0546596, + 0.0634012, + 0.0613207, + 0.0653995, + 0.0450429, + 0.065412, + 0.0585679, + 0.0530456, + 0.0484401, + 0.0694564, + 0.0684182, + 0.117808, + 0.126668, + 0.142361, + 0.154379, + 0.102089, + 0.114404, + 0.160111, + 0.178321, + 0.0964375, + 0.131446, + 0.115852, + 0.326339 + ) + +TrgLayer_weights = cms.vdouble(0.0, + 0.0183664, + 0., + 0.0305622, + 0., + 0.0162589, + 0., + 0.0143918, + 0., + 0.0134475, + 0., + 0.0185754, + 0., + 0.0204934, + 0., + 0.016901, + 0., + 0.0207958, + 0., + 0.0167985, + 0., + 0.0238385, + 0., + 0.0301146, + 0., + 0.0274622, + 0., + 0.0468671, + 0., + 0.078819, + 0.0555831, + 0.0609312, + 0.0610768, + 0.0657626, + 0.0465653, + 0.0629383, + 0.0610061, + 0.0517326, + 0.0492882, + 0.0699336, + 0.0673457, + 0.119896, + 0.125327, + 0.143235, + 0.153295, + 0.104777, + 0.109345, + 0.161386, + 0.174656, + 0.108053, + 0.121674, + 0.1171, + 0.328053 + ) + +TrgLayer_dEdX_weights = cms.vdouble(0.0, # there is no layer zero + 8.603+8.0675, # Mev + 0., + 8.0675*2, + 0., + 8.0675*2, + 0., + 8.0675*2., + 0., + 8.0675+8.9515, + 0., + 10.135*2, + 0., + 10.135*2., + 0., + 10.135*2., + 0., + 10.135*2., + 0., + 10.135+11.682, + 0., + 13.654*2., + 0., + 13.654*2., + 0., + 13.654*2., + 0., + 13.654+38.2005, + 0., + 55.0265, + 49.871, + 49.871, + 49.871, + 49.871, + 49.871, + 49.871, + 49.871, + 49.871, + 49.871, + 49.871, + 62.005, + 83.1675, + 92.196, + 92.196, + 92.196, + 92.196, + 92.196, + 92.196, + 92.196, + 92.196, + 92.196, + 92.196, + 46.098) diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerGeometryESProducer_cfi.py b/L1Trigger/L1THGCal/python/hgcalTriggerGeometryESProducer_cfi.py index 57f112097ff0b..b1b1e4ee7eff7 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerGeometryESProducer_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalTriggerGeometryESProducer_cfi.py @@ -1,13 +1,31 @@ import FWCore.ParameterSet.Config as cms +disconnectedTriggerLayers = [ + 2, + 4, + 6, + 8, + 10, + 12, + 14, + 16, + 18, + 20, + 22, + 24, + 26, + 28 + ] + geometry = cms.PSet( TriggerGeometryName = cms.string('HGCalTriggerGeometryHexLayerBasedImp1'), L1TCellsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_8inch_aligned_192_432_V8_0.txt"), - L1TModulesMapping = cms.FileInPath("L1Trigger/L1THGCal/data/panel_mapping_60deg_6mod_0.txt"), + L1TModulesMapping = cms.FileInPath("L1Trigger/L1THGCal/data/panel_mapping_tdr_0.txt"), L1TCellNeighborsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_8inch_aligned_192_432_0.txt"), L1TCellsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_BH_3x3_30deg_0.txt"), L1TCellNeighborsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_BH_3x3_30deg_0.txt"), - DisconnectedModules = cms.vuint32(0) + DisconnectedModules = cms.vuint32(0), + DisconnectedLayers = cms.vuint32(disconnectedTriggerLayers) ) hgcalTriggerGeometryESProducer = cms.ESProducer( diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py index 9d2cae41a06f8..2f6f2e5b818a1 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalTriggerNtuples_cfi.py @@ -1,19 +1,37 @@ import FWCore.ParameterSet.Config as cms +import SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi as digiparam +import RecoLocalCalo.HGCalRecProducers.HGCalUncalibRecHit_cfi as recoparam +import RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi as recocalibparam +import hgcalLayersCalibrationCoefficients_cfi as layercalibparam + + +fcPerMip = recoparam.HGCalUncalibRecHit.HGCEEConfig.fCPerMIP +keV2fC = digiparam.hgceeDigitizer.digiCfg.keV2fC +layerWeights = layercalibparam.TrgLayer_dEdX_weights +thicknessCorrections = recocalibparam.HGCalRecHit.thicknessCorrection ntuple_event = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleEvent') ) + +from FastSimulation.Event.ParticleFilter_cfi import ParticleFilterBlock +PartFilterConfig = ParticleFilterBlock.ParticleFilter.copy() +PartFilterConfig.protonEMin = cms.double(100000) +PartFilterConfig.etaMax = cms.double(3.1) + ntuple_gen = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleGen'), - GenParticles = cms.InputTag('genParticles') + GenParticles = cms.InputTag('genParticles'), + GenPU = cms.InputTag('addPileupInfo'), + particleFilter = PartFilterConfig ) ntuple_gentau = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleGenTau'), GenParticles = cms.InputTag('genParticles'), - isPythia8 = cms.bool(True) + isPythia8 = cms.bool(False) ) ntuple_genjet = cms.PSet( @@ -34,12 +52,24 @@ ntuple_triggercells = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCTriggerCells'), - TriggerCells = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:calibratedTriggerCells') + TriggerCells = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:calibratedTriggerCells'), + Multiclusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster3D'), + eeSimHits = cms.InputTag('g4SimHits:HGCHitsEE'), + fhSimHits = cms.InputTag('g4SimHits:HGCHitsHEfront'), + bhSimHits = cms.InputTag('g4SimHits:HcalHits'), + FillSimEnergy = cms.bool(False), + fcPerMip = fcPerMip, + keV2fC = keV2fC, + layerWeights = layerWeights, + thicknessCorrections = thicknessCorrections, + FilterCellsInMulticlusters = cms.bool(True) ) ntuple_clusters = cms.PSet( NtupleName = cms.string('HGCalTriggerNtupleHGCClusters'), - Clusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster2D') + Clusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster2D'), + Multiclusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster3D'), + FilterClustersInMulticlusters = cms.bool(True) ) ntuple_multicluster = cms.PSet( @@ -47,12 +77,18 @@ Multiclusters = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:cluster3D') ) +ntuple_panels = cms.PSet( + NtupleName = cms.string('HGCalTriggerNtupleHGCPanels'), + TriggerCells = cms.InputTag('hgcalTriggerPrimitiveDigiProducer:calibratedTriggerCells') +) + hgcalTriggerNtuplizer = cms.EDAnalyzer( "HGCalTriggerNtupleManager", Ntuples = cms.VPSet( ntuple_event, ntuple_gen, ntuple_genjet, + ntuple_gentau, ntuple_digis, ntuple_triggercells, ntuple_clusters, diff --git a/L1Trigger/L1THGCal/python/hgcalTriggerPrimitiveDigiProducer_cfi.py b/L1Trigger/L1THGCal/python/hgcalTriggerPrimitiveDigiProducer_cfi.py index 24395ff0c5787..603feb4bb513e 100644 --- a/L1Trigger/L1THGCal/python/hgcalTriggerPrimitiveDigiProducer_cfi.py +++ b/L1Trigger/L1THGCal/python/hgcalTriggerPrimitiveDigiProducer_cfi.py @@ -2,6 +2,7 @@ import SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi as digiparam import RecoLocalCalo.HGCalRecProducers.HGCalUncalibRecHit_cfi as recoparam import RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi as recocalibparam +import hgcalLayersCalibrationCoefficients_cfi as layercalibparam # Digitization parameters adcSaturation_fC = digiparam.hgceeDigitizer.digiCfg.feCfg.adcSaturation_fC @@ -14,7 +15,7 @@ # Reco calibration parameters fCPerMIPee = recoparam.HGCalUncalibRecHit.HGCEEConfig.fCPerMIP fCPerMIPfh = recoparam.HGCalUncalibRecHit.HGCHEFConfig.fCPerMIP -layerWeights = recocalibparam.HGCalRecHit.layerWeights +layerWeights = layercalibparam.TrgLayer_dEdX_weights thicknessCorrection = recocalibparam.HGCalRecHit.thicknessCorrection # Parameters used in several places @@ -31,13 +32,13 @@ fe_codec = cms.PSet( CodecName = cms.string('HGCalTriggerCellThresholdCodec'), CodecIndex = cms.uint32(2), MaxCellsInModule = cms.uint32(288), - DataLength = cms.uint32(16), + DataLength = cms.uint32(20), linLSB = cms.double(triggerCellLsbBeforeCompression), linnBits = cms.uint32(16), triggerCellTruncationBits = cms.uint32(triggerCellTruncationBits), NData = cms.uint32(999), - TCThreshold_fC = cms.double(1.), - TCThresholdBH_MIP = cms.double(1.), + TCThreshold_fC = cms.double(0.), + TCThresholdBH_MIP = cms.double(0.), #take the following parameters from the digitization config file adcsaturation = adcSaturation_fC, adcnBits = adcNbits, @@ -55,21 +56,32 @@ dEdXweights = layerWeights, thickCorr = cms.double(thicknessCorrection_200) ) + C2d_parValues = cms.PSet( seeding_threshold_silicon = cms.double(5), # MipT seeding_threshold_scintillator = cms.double(5), # MipT clustering_threshold_silicon = cms.double(2), # MipT clustering_threshold_scintillator = cms.double(2), # MipT - dR_cluster = cms.double(3.), # in cm - clusterType = cms.string('NNC2d') # clustering type: dRC2d--> Geometric-dR clustering; NNC2d-->Nearest Neighbors clustering + clusterType = cms.string('NNC2d'), + applyLayerCalibration = cms.bool(True), + layerWeights = layercalibparam.TrgLayer_weights, + # Parameters not used by this clustering + dR_cluster=cms.double(0.), + calibSF_cluster=cms.double(0.) ) C3d_parValues = cms.PSet( dR_multicluster = cms.double(0.01), # dR in normalized plane used to clusterize C2d minPt_multicluster = cms.double(0.5), # minimum pt of the multicluster (GeV) - calibSF_multicluster = cms.double(1.084) - ) + type_multicluster = cms.string('dRC3d'), + # Parameters not used by this clustering + dist_dbscan_multicluster=cms.double(0.), + minN_dbscan_multicluster=cms.uint32(0) +) + cluster_algo = cms.PSet( AlgorithmName = cms.string('HGCClusterAlgoThreshold'), FECodec = fe_codec.clone(), calib_parameters = calib_parValues.clone(), + triggercell_threshold_silicon = cms.double(2.), # MipT + triggercell_threshold_scintillator = cms.double(2.), # MipT C2d_parameters = C2d_parValues.clone(), C3d_parameters = C3d_parValues.clone() ) diff --git a/L1Trigger/L1THGCal/src/HGCalTriggerGeometryGenericMapping.cc b/L1Trigger/L1THGCal/src/HGCalTriggerGeometryGenericMapping.cc index dc0bf68bffff8..138b06f9ddd43 100644 --- a/L1Trigger/L1THGCal/src/HGCalTriggerGeometryGenericMapping.cc +++ b/L1Trigger/L1THGCal/src/HGCalTriggerGeometryGenericMapping.cc @@ -1,5 +1,6 @@ #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryGenericMapping.h" + using namespace HGCalTriggerGeometry; namespace { @@ -131,3 +132,10 @@ disconnectedModule(const unsigned module_id) const { return false; } + +unsigned +HGCalTriggerGeometryGenericMapping:: +triggerLayer(const unsigned id) const { + return HGCalDetId(id).layer(); +} + diff --git a/L1Trigger/L1THGCal/src/HGCalTriggerTools.cc b/L1Trigger/L1THGCal/src/HGCalTriggerTools.cc new file mode 100644 index 0000000000000..0df4ea772567c --- /dev/null +++ b/L1Trigger/L1THGCal/src/HGCalTriggerTools.cc @@ -0,0 +1,158 @@ +#include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h" + + +#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h" + +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "DataFormats/HcalDetId/interface/HcalDetId.h" +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" +#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" + +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/Framework/interface/EventSetup.h" + + +namespace { + constexpr char hgcalee_sens[] = "HGCalEESensitive"; + constexpr char hgcalfh_sens[] = "HGCalHESiliconSensitive"; + + constexpr std::float_t idx_to_thickness = std::float_t(100.0); + + template + inline void check_ddd(const DDD* ddd) { + if( nullptr == ddd ) { + throw cms::Exception("hgcal::HGCalTriggerTools") + << "DDDConstants not accessibl to hgcal::HGCalTriggerTools!"; + } + } + + template + inline void check_geom(const GEOM* geom) { + if( nullptr == geom ) { + throw cms::Exception("hgcal::HGCalTriggerTools") + << "Geometry not provided yet to hgcal::HGCalTriggerTools!"; + } + } + + inline const HcalDDDRecConstants* get_ddd(const CaloSubdetectorGeometry* geom, + const HcalDetId& detid) { + const HcalGeometry* hc = static_cast(geom); + const HcalDDDRecConstants* ddd = hc->topology().dddConstants(); + check_ddd(ddd); + return ddd; + } + + inline const HGCalDDDConstants* get_ddd(const CaloSubdetectorGeometry* geom, + const HGCalDetId& detid) { + const HGCalGeometry* hg = static_cast(geom); + const HGCalDDDConstants* ddd = &(hg->topology().dddConstants()); + check_ddd(ddd); + return ddd; + } + +} + +void HGCalTriggerTools::setEventSetup(const edm::EventSetup& es) { + + edm::ESHandle triggerGeometry_; + es.get().get(triggerGeometry_); + geom_ = triggerGeometry_.product(); + fhOffset_ = (geom_->eeTopology().dddConstants()).layers(true); + bhOffset_ = fhOffset_ + (geom_->fhTopology().dddConstants()).layers(true); +} + +GlobalPoint HGCalTriggerTools::getTCPosition(const DetId& id) const { + if(id.det() == DetId::Hcal) { + throw cms::Exception("hgcal::HGCalTriggerTools") + << "method getTCPosition called for DetId not belonging to a TC"; + // FIXME: this would actually need a better test...but at the moment I can not think to anything better + // to distinguish a TC detId + } + + GlobalPoint position = geom_->getTriggerCellPosition(id); + return position; +} + + +unsigned int HGCalTriggerTools::getLayer(const DetId& id) const { + unsigned int layer = std::numeric_limits::max(); + if( id.det() == DetId::Forward) { + const HGCalDetId hid(id); + layer = hid.layer(); + } else if( id.det() == DetId::Hcal && id.subdetId() == HcalEndcap) { + const HcalDetId hcid(id); + layer = hcid.depth(); + } + return layer; +} + +unsigned int HGCalTriggerTools::getLayerWithOffset(const DetId& id) const { + unsigned int layer = getLayer(id); + if( id.det() == DetId::Forward && id.subdetId() == HGCHEF ) { + layer += fhOffset_; + } else if( (id.det() == DetId::Hcal && id.subdetId() == HcalEndcap) || + (id.det() == DetId::Forward && id.subdetId() == HGCHEB) ) { + layer += bhOffset_; + } + return layer; +} + +float HGCalTriggerTools::getEta(const GlobalPoint& position, const float& vertex_z) const { + GlobalPoint corrected_position = GlobalPoint(position.x(), position.y(), position.z()-vertex_z); + return corrected_position.eta(); +} + +float HGCalTriggerTools::getTCEta(const DetId& id, const float& vertex_z) const { + GlobalPoint position = getTCPosition(id); + return getEta(position, vertex_z); +} + +float HGCalTriggerTools::getPhi(const GlobalPoint& position) const { + float phi = atan2(position.y(),position.x()); + return phi; +} + +float HGCalTriggerTools::getTCPhi(const DetId& id) const { + GlobalPoint position = getTCPosition(id); + return getPhi(position); +} + +float HGCalTriggerTools::getPt(const GlobalPoint& position, const float& hitEnergy, const float& vertex_z) const { + float eta = getEta(position, vertex_z); + float pt = hitEnergy / cosh(eta); + return pt; +} + +float HGCalTriggerTools::getTCPt(const DetId& id, const float& hitEnergy, const float& vertex_z) const { + GlobalPoint position = getTCPosition(id); + return getPt(position, hitEnergy, vertex_z); +} + +float HGCalTriggerTools::getLayerZ(const unsigned& layerWithOffset) const { + int subdet = ForwardSubdetector::HGCEE; + unsigned offset = 0; + if(layerWithOffset > lastLayerEE() && layerWithOffset <= lastLayerFH()) { + subdet = ForwardSubdetector::HGCHEF; + offset = lastLayerEE(); + } else if(layerWithOffset > lastLayerFH()) { + subdet = HcalSubdetector::HcalEndcap; + offset = lastLayerFH(); + } + return getLayerZ(subdet, layerWithOffset - offset); +} + +float HGCalTriggerTools::getLayerZ(const int& subdet, const unsigned& layer) const { + const unsigned heOffset = 7; + float layerGlobalZ = 0.; + if(subdet == ForwardSubdetector::HGCEE) { + layerGlobalZ = geom_->eeTopology().dddConstants().waferZ(layer, true); + } else if(subdet == ForwardSubdetector::HGCHEF) { + layerGlobalZ = geom_->fhTopology().dddConstants().waferZ(layer, true); + } else if(subdet == HcalSubdetector::HcalEndcap || subdet == ForwardSubdetector::HGCHEB) { + layerGlobalZ = geom_->bhTopology().dddConstants()->getRZ(HcalSubdetector::HcalEndcap, layer+heOffset); + } + return layerGlobalZ; +} diff --git a/L1Trigger/L1THGCal/src/be_algorithms/HGCalClusteringImpl.cc b/L1Trigger/L1THGCal/src/be_algorithms/HGCalClusteringImpl.cc index d64b5e859f4bc..f8174971f9e92 100644 --- a/L1Trigger/L1THGCal/src/be_algorithms/HGCalClusteringImpl.cc +++ b/L1Trigger/L1THGCal/src/be_algorithms/HGCalClusteringImpl.cc @@ -11,14 +11,17 @@ HGCalClusteringImpl::HGCalClusteringImpl(const edm::ParameterSet & conf): scintillatorSeedThreshold_(conf.getParameter("seeding_threshold_scintillator")), scintillatorTriggerCellThreshold_(conf.getParameter("clustering_threshold_scintillator")), dr_(conf.getParameter("dR_cluster")), - clusteringAlgorithmType_(conf.getParameter("clusterType")) + clusteringAlgorithmType_(conf.getParameter("clusterType")), + calibSF_(conf.getParameter("calibSF_cluster")), + layerWeights_(conf.getParameter< std::vector >("layerWeights")), + applyLayerWeights_(conf.getParameter< bool >("applyLayerCalibration")) { edm::LogInfo("HGCalClusterParameters") << "C2d Clustering Algorithm selected : " << clusteringAlgorithmType_ ; edm::LogInfo("HGCalClusterParameters") << "C2d silicon seeding Thr: " << siliconSeedThreshold_ ; edm::LogInfo("HGCalClusterParameters") << "C2d silicon clustering Thr: " << siliconTriggerCellThreshold_ ; edm::LogInfo("HGCalClusterParameters") << "C2d scintillator seeding Thr: " << scintillatorSeedThreshold_ ; edm::LogInfo("HGCalClusterParameters") << "C2d scintillator clustering Thr: " << scintillatorTriggerCellThreshold_ ; - + edm::LogInfo("HGCalClusterParameters") << "C2d global calibration factor: " << calibSF_; } @@ -43,15 +46,15 @@ bool HGCalClusteringImpl::isPertinent( const l1t::HGCalTriggerCell & tc, } -void HGCalClusteringImpl::clusterizeDR( const edm::PtrVector & triggerCellsPtrs, - l1t::HGCalClusterBxCollection & clusters +void HGCalClusteringImpl::clusterizeDR( const std::vector> & triggerCellsPtrs, + l1t::HGCalClusterBxCollection & clusters ){ bool isSeed[triggerCellsPtrs.size()]; /* search for cluster seeds */ int itc(0); - for( edm::PtrVector::const_iterator tc = triggerCellsPtrs.begin(); tc != triggerCellsPtrs.end(); ++tc,++itc ){ + for( std::vector>::const_iterator tc = triggerCellsPtrs.begin(); tc != triggerCellsPtrs.end(); ++tc,++itc ){ double seedThreshold = ((*tc)->subdetId()==HGCHEB ? scintillatorSeedThreshold_ : siliconSeedThreshold_); isSeed[itc] = ( (*tc)->mipPt() > seedThreshold) ? true : false; } @@ -60,13 +63,17 @@ void HGCalClusteringImpl::clusterizeDR( const edm::PtrVector clustersTmp; itc=0; - for( edm::PtrVector::const_iterator tc = triggerCellsPtrs.begin(); tc != triggerCellsPtrs.end(); ++tc,++itc ){ + for( std::vector>::const_iterator tc = triggerCellsPtrs.begin(); tc != triggerCellsPtrs.end(); ++tc,++itc ){ double threshold = ((*tc)->subdetId()==HGCHEB ? scintillatorTriggerCellThreshold_ : siliconTriggerCellThreshold_); if( (*tc)->mipPt() < threshold ){ continue; } - /* searching for TC near the center of the cluster */ + /* + searching for TC near the center of the cluster + ToBeFixed: if a tc is not a seed, but could be potencially be part of a cluster generated by a late seed, + the tc will not be clusterized + */ int iclu=0; vector tcPertinentClusters; for( const auto& clu : clustersTmp){ @@ -99,9 +106,12 @@ void HGCalClusteringImpl::clusterizeDR( const edm::PtrVector & triggerCellsPtrs, +void HGCalClusteringImpl::triggerCellReshuffling( const std::vector> & triggerCellsPtrs, std::array< std::array>, kLayers_>,kNSides_> & reshuffledTriggerCells ){ @@ -141,9 +151,9 @@ void HGCalClusteringImpl::mergeClusters( l1t::HGCalCluster & main_cluster, const l1t::HGCalCluster & secondary_cluster ) const { - const edm::PtrVector& pertinentTC = secondary_cluster.constituents(); + const std::vector>& pertinentTC = secondary_cluster.constituents(); - for( edm::PtrVector::iterator tc = pertinentTC.begin(); tc != pertinentTC.end(); ++tc ){ + for( std::vector>::const_iterator tc = pertinentTC.begin(); tc != pertinentTC.end(); ++tc ){ main_cluster.addConstituent(*tc); } @@ -248,15 +258,16 @@ void HGCalClusteringImpl::NNKernel( const std::vector & triggerCellsPtrs, - l1t::HGCalClusterBxCollection & clusters, - const HGCalTriggerGeometryBase & triggerGeometry +void HGCalClusteringImpl::clusterizeNN( const std::vector> & triggerCellsPtrs, + l1t::HGCalClusterBxCollection & clusters, + const HGCalTriggerGeometryBase & triggerGeometry ){ std::array< std::array< std::vector >,kLayers_>,kNSides_> reshuffledTriggerCells; @@ -270,3 +281,201 @@ void HGCalClusteringImpl::clusterizeNN( const edm::PtrVector> & triggerCellsPtrs, + l1t::HGCalClusterBxCollection & clusters, + const HGCalTriggerGeometryBase & triggerGeometry + ){ + + bool isSeed[triggerCellsPtrs.size()]; + std::vector seedPositions; + seedPositions.reserve( triggerCellsPtrs.size() ); + + /* search for cluster seeds */ + int itc(0); + for( std::vector>::const_iterator tc = triggerCellsPtrs.begin(); tc != triggerCellsPtrs.end(); ++tc,++itc ){ + + double seedThreshold = ((*tc)->subdetId()==HGCHEB ? scintillatorSeedThreshold_ : siliconSeedThreshold_); + + /* decide if is a seed, if yes store the position into of triggerCellsPtrs */ + isSeed[itc] = ( (*tc)->mipPt() > seedThreshold) ? true : false; + if( isSeed[itc] ) { + + seedPositions.push_back( itc ); + + /* remove tc from the seed vector if is a NN of an other seed*/ + for( auto pos : seedPositions ){ + if( ( (*tc)->position() - triggerCellsPtrs[pos]->position() ).mag() < dr_ ){ + if( this->areTCneighbour( (*tc)->detId(), triggerCellsPtrs[pos]->detId(), triggerGeometry ) ) + { + isSeed[itc] = false; + seedPositions.pop_back(); + } + } + } + } + + } + + /* clustering the TCs */ + std::vector clustersTmp; + + // every seed generates a cluster + for( auto pos : seedPositions ) { + clustersTmp.emplace_back( triggerCellsPtrs[pos] ); + } + + /* add the tc to the clusters */ + itc=0; + for( std::vector>::const_iterator tc = triggerCellsPtrs.begin(); tc != triggerCellsPtrs.end(); ++tc,++itc ){ + + /* get the correct threshold for the different part of the detector */ + double threshold = ((*tc)->subdetId()==HGCHEB ? scintillatorTriggerCellThreshold_ : siliconTriggerCellThreshold_); + + /* continue if not passing the threshold */ + if( (*tc)->mipPt() < threshold ) continue; + if( isSeed[itc] ) continue; //No sharing of seeds between clusters (TBC) + + /* searching for TC near the centre of the cluster */ + std::vector tcPertinentClusters; + unsigned iclu(0); + + for ( const auto& clu : clustersTmp ) { + if( this->isPertinent(**tc, clu, dr_) ){ + tcPertinentClusters.push_back( iclu ); + } + ++iclu; + } + + if ( tcPertinentClusters.size() == 0 ) { + continue; + } + else if( tcPertinentClusters.size() == 1 ) { + clustersTmp.at( tcPertinentClusters.at(0) ).addConstituent( *tc ); + } + else { + + /* calculate the fractions */ + double totMipt = 0; + for( auto clu : tcPertinentClusters ){ + totMipt += clustersTmp.at( clu ).seedMipPt(); + } + + for( auto clu : tcPertinentClusters ){ + double seedMipt = clustersTmp.at( clu ).seedMipPt(); + clustersTmp.at( clu ).addConstituent( *tc, true, seedMipt/totMipt ); + } + } + } + + /* store clusters in the persistent collection */ + clusters.resize(0, clustersTmp.size()); + for( unsigned i(0); iremoveUnconnectedTCinCluster( clustersTmp.at(i), triggerGeometry ); + calibratePt( clustersTmp.at(i) ); + clusters.set( 0, i, clustersTmp.at(i) ); + } + +} + + +bool HGCalClusteringImpl::areTCneighbour(uint32_t detIDa, uint32_t detIDb, const HGCalTriggerGeometryBase & triggerGeometry + ){ + + const auto neighbors = triggerGeometry.getNeighborsFromTriggerCell( detIDa ); + + if( neighbors.find( detIDb ) != neighbors.end() ) return true; + + return false; + +} + + +void HGCalClusteringImpl::removeUnconnectedTCinCluster( l1t::HGCalCluster & cluster, const HGCalTriggerGeometryBase & triggerGeometry ) { + + /* get the constituents and the centre of the seed tc (considered as the first of the constituents) */ + const std::vector>& constituents = cluster.constituents(); + Basic3DVector seedCentre( constituents[0]->position() ); + + /* distances from the seed */ + vector,float>> distances; + for( const auto & tc : constituents ) + { + Basic3DVector tcCentre( tc->position() ); + float distance = ( seedCentre - tcCentre ).mag(); + distances.push_back( pair,float>( tc, distance ) ); + } + + /* sorting (needed in order to be sure that we are skipping any tc) */ + /* FIXME: better sorting needed!!! */ + std::sort( distances.begin(), distances.end(), distanceSorter ); + + + /* checking if the tc is connected to the seed */ + bool toRemove[constituents.size()]; + toRemove[0] = false; // this is the seed + for( unsigned itc=1; itc& tcToStudy = distances[itc].first; + + /* compare with the tc in the cluster */ + for( unsigned itc_ref=0; itc_refdetId(), distances.at( itc_ref ).first->detId(), triggerGeometry ) ) { + toRemove[itc] = false; + break; + } + } + } + + } + + + /* remove the unconnected TCs */ + for( unsigned i=0; i("dR_multicluster")), ptC3dThreshold_(conf.getParameter("minPt_multicluster")), - calibSF_(conf.getParameter("calibSF_multicluster")) + multiclusterAlgoType_(conf.getParameter("type_multicluster")), + distDbscan_(conf.getParameter("dist_dbscan_multicluster")), + minNDbscan_(conf.getParameter("minN_dbscan_multicluster")) { edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster dR for Near Neighbour search: " << dr_; edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster minimum transverse-momentum: " << ptC3dThreshold_; - edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster global calibration factor: " << calibSF_; - + edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster DBSCAN Clustering distance: " << distDbscan_; + edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster clustering min number of subclusters: " << minNDbscan_; + edm::LogInfo("HGCalMulticlusterParameters") << "Multicluster type of multiclustering algortihm: " << multiclusterAlgoType_; } @@ -32,14 +38,47 @@ bool HGCalMulticlusteringImpl::isPertinent( const l1t::HGCalCluster & clu, } -void HGCalMulticlusteringImpl::clusterize( const edm::PtrVector & clustersPtrs, - l1t::HGCalMulticlusterBxCollection & multiclusters) +void HGCalMulticlusteringImpl::findNeighbor( const std::vector>& rankedList, + unsigned int searchInd, + const std::vector>& clustersPtrs, + std::vector& neighbors + ){ + + if(clustersPtrs.size() <= searchInd || clustersPtrs.size() < rankedList.size()){ + throw cms::Exception("IndexOutOfBound: clustersPtrs in 'findNeighbor'"); + } + + for(unsigned int ind = searchInd+1; ind < rankedList.size() && fabs(rankedList.at(ind).second - rankedList.at(searchInd).second) < distDbscan_ ; ind++){ + + if(clustersPtrs.size() <= rankedList.at(ind).first){ + throw cms::Exception("IndexOutOfBound: clustersPtrs in 'findNeighbor'"); + + } else if(((*(clustersPtrs[rankedList.at(ind).first])).centreProj() - (*(clustersPtrs[rankedList.at(searchInd).first])).centreProj()).mag() < distDbscan_){ + neighbors.push_back(ind); + } + } + + for(unsigned int ind = 0; ind < searchInd && fabs(rankedList.at(searchInd).second - rankedList.at(ind).second) < distDbscan_ ; ind++){ + + if(clustersPtrs.size() <= rankedList.at(ind).first){ + throw cms::Exception("IndexOutOfBound: clustersPtrs in 'findNeighbor'"); + + } else if(((*(clustersPtrs[rankedList.at(ind).first])).centreProj() - (*(clustersPtrs[rankedList.at(searchInd).first])).centreProj()).mag() < distDbscan_){ + neighbors.push_back(ind); + } + } +} + + +void HGCalMulticlusteringImpl::clusterizeDR( const std::vector> & clustersPtrs, + l1t::HGCalMulticlusterBxCollection & multiclusters, + const HGCalTriggerGeometryBase & triggerGeometry) { - + std::vector multiclustersTmp; int iclu = 0; - for(edm::PtrVector::const_iterator clu = clustersPtrs.begin(); clu != clustersPtrs.end(); ++clu, ++iclu){ + for(std::vector>::const_iterator clu = clustersPtrs.begin(); clu != clustersPtrs.end(); ++clu, ++iclu){ int imclu=0; vector tcPertinentMulticlusters; @@ -70,15 +109,138 @@ void HGCalMulticlusteringImpl::clusterize( const edm::PtrVector> &pertinentClu = multiclustersTmp.at(i).constituents(); + for( std::vector>::const_iterator it_clu=pertinentClu.begin(); it_clupt(); + math::PtEtaPhiMLorentzVector multiclusterP4( sumPt, + multiclustersTmp.at(i).centre().eta(), + multiclustersTmp.at(i).centre().phi(), + 0. ); + multiclustersTmp.at(i).setP4( multiclusterP4 ); + + if( multiclustersTmp.at(i).pt() > ptC3dThreshold_ ){ + + //compute shower shape + multiclustersTmp.at(i).set_showerLength(shape_.showerLength(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_coreShowerLength(shape_.coreShowerLength(multiclustersTmp.at(i), triggerGeometry)); + multiclustersTmp.at(i).set_firstLayer(shape_.firstLayer(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_maxLayer(shape_.maxLayer(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaEtaEtaTot(shape_.sigmaEtaEtaTot(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaEtaEtaMax(shape_.sigmaEtaEtaMax(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaPhiPhiTot(shape_.sigmaPhiPhiTot(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaPhiPhiMax(shape_.sigmaPhiPhiMax(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaZZ(shape_.sigmaZZ(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaRRTot(shape_.sigmaRRTot(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaRRMax(shape_.sigmaRRMax(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaRRMean(shape_.sigmaRRMean(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_eMax(shape_.eMax(multiclustersTmp.at(i))); + multiclusters.push_back( 0, multiclustersTmp.at(i)); } } } +void HGCalMulticlusteringImpl::clusterizeDBSCAN( const std::vector> & clustersPtrs, + l1t::HGCalMulticlusterBxCollection & multiclusters, + const HGCalTriggerGeometryBase & triggerGeometry) +{ + + std::vector multiclustersTmp; + l1t::HGCalMulticluster mcluTmp; + std::vector visited(clustersPtrs.size(),false); + std::vector merged (clustersPtrs.size(),false); + std::vector> rankedList; + rankedList.reserve(clustersPtrs.size()); + std::vector> neighborList; + neighborList.reserve(clustersPtrs.size()); + + int iclu = 0, imclu = 0, neighNo; + double dist = 0.; + + for(std::vector>::const_iterator clu = clustersPtrs.begin(); clu != clustersPtrs.end(); ++clu, ++iclu){ + dist = (*clu)->centreProj().mag()*HGCalDetId((*clu)->detId()).zside(); + rankedList.push_back(std::make_pair(iclu,dist)); + } + iclu = 0; + std::sort(rankedList.begin(), rankedList.end(), [](auto &left, auto &right) { + return left.second < right.second; + }); + + for(auto cluRanked: rankedList){ + std::vector neighbors; + + if(!visited.at(iclu)){ + visited.at(iclu) = true; + findNeighbor(rankedList, iclu, clustersPtrs, neighbors); + neighborList.push_back(std::move(neighbors)); + + if(neighborList.at(iclu).size() >= minNDbscan_) { + multiclustersTmp.emplace_back( clustersPtrs[cluRanked.first] ); + merged.at(iclu) = true; + /* dynamic range loop: range-based loop syntax cannot be employed */ + for(unsigned int neighInd = 0; neighInd < neighborList.at(iclu).size(); neighInd++){ + neighNo = neighborList.at(iclu).at(neighInd); + + if(!visited.at(neighNo)){ + visited.at(neighNo) = true; + std::vector secNeighbors; + findNeighbor(rankedList, neighNo,clustersPtrs, secNeighbors); + multiclustersTmp.at(imclu).addConstituent( clustersPtrs[rankedList.at(neighNo).first]); + merged.at(neighNo) = true; + + if(secNeighbors.size() >= minNDbscan_){ + neighborList.at(iclu).insert(neighborList.at(iclu).end(), secNeighbors.begin(), secNeighbors.end()); + } + + } else if(!merged.at(neighNo) ){ + merged.at(neighNo) = true; + multiclustersTmp.at(imclu).addConstituent( clustersPtrs[rankedList.at(neighNo).first] ); + } + } + imclu++; + } + } + else neighborList.push_back(std::move(neighbors)); + iclu++; + } + /* making the collection of multiclusters */ + for( unsigned i(0); i> &pertinentClu = multiclustersTmp.at(i).constituents(); + for( std::vector>::const_iterator it_clu=pertinentClu.begin(); it_clupt(); + + math::PtEtaPhiMLorentzVector multiclusterP4( sumPt, + multiclustersTmp.at(i).centre().eta(), + multiclustersTmp.at(i).centre().phi(), + 0. ); + multiclustersTmp.at(i).setP4( multiclusterP4 ); + + + if( multiclustersTmp.at(i).pt() > ptC3dThreshold_ ){ + + //compute shower shape + multiclustersTmp.at(i).set_showerLength(shape_.showerLength(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_firstLayer(shape_.firstLayer(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaEtaEtaTot(shape_.sigmaEtaEtaTot(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaEtaEtaMax(shape_.sigmaEtaEtaMax(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaPhiPhiTot(shape_.sigmaPhiPhiTot(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaPhiPhiMax(shape_.sigmaPhiPhiMax(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaZZ(shape_.sigmaZZ(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaRRTot(shape_.sigmaRRTot(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_sigmaRRMax(shape_.sigmaRRMax(multiclustersTmp.at(i))); + multiclustersTmp.at(i).set_eMax(shape_.eMax(multiclustersTmp.at(i))); + + multiclusters.push_back( 0, multiclustersTmp.at(i)); + } + } + +} diff --git a/L1Trigger/L1THGCal/src/be_algorithms/HGCalShowerShape.cc b/L1Trigger/L1THGCal/src/be_algorithms/HGCalShowerShape.cc new file mode 100644 index 0000000000000..9cc6ed9c3af30 --- /dev/null +++ b/L1Trigger/L1THGCal/src/be_algorithms/HGCalShowerShape.cc @@ -0,0 +1,549 @@ +#include "L1Trigger/L1THGCal/interface/be_algorithms/HGCalShowerShape.h" +#include "TMath.h" +#include +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" +#include "DataFormats/L1THGCal/interface/HGCalCluster.h" +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +#include +#include +#include +#include +#include + + +int HGCalShowerShape::HGC_layer(const uint32_t subdet, const uint32_t layer) const { + + int hgclayer = -1; + if(subdet==HGCEE) hgclayer=layer;//EE + else if(subdet==HGCHEF) hgclayer=layer+kLayersEE_;//FH + else if(subdet==HGCHEB) hgclayer=layer+kLayersEE_+kLayersFH_;//BH + + return hgclayer; + +} + +//Compute energy-weighted mean of any variable X in the cluster + +float HGCalShowerShape::meanX(const std::vector >& energy_X_tc) const { + + float Etot = 0; + float X_sum = 0; + + for(const auto& energy_X : energy_X_tc){ + + X_sum += energy_X.first*energy_X.second; + Etot += energy_X.first; + + } + + float X_mean = 0; + if(Etot>0) X_mean = X_sum/Etot; + return X_mean; + +} + +//Compute energy-weighted RMS of any variable X in the cluster + +float HGCalShowerShape::sigmaXX(const std::vector >& energy_X_tc, const float X_cluster) const { + + float Etot = 0; + float deltaX2_sum = 0; + + for(const auto& energy_X : energy_X_tc){ + + deltaX2_sum += energy_X.first * pow(energy_X.second - X_cluster,2); + Etot += energy_X.first; + + } + + float X_MSE = 0; + if (Etot>0) X_MSE=deltaX2_sum/Etot; + float X_RMS=sqrt(X_MSE); + return X_RMS; + +} + + +//Compute energy-weighted RMS of any variable X in the cluster +//Extra care needed because of deltaPhi + +float HGCalShowerShape::sigmaPhiPhi(const std::vector >& energy_phi_tc, const float phi_cluster) const { + + float Etot = 0; + float deltaphi2_sum = 0; + + for(const auto& energy_phi : energy_phi_tc){ + + deltaphi2_sum += energy_phi.first * pow(deltaPhi(energy_phi.second,phi_cluster),2); + Etot += energy_phi.first; + + } + + float phi_MSE = 0; + if (Etot>0) phi_MSE=deltaphi2_sum/Etot; + float Spp=sqrt(phi_MSE); + return Spp; + +} + + + +int HGCalShowerShape::firstLayer(const l1t::HGCalMulticluster& c3d) const { + + const std::vector>& clustersPtrs = c3d.constituents(); + + int firstLayer=999; + + for(const auto& clu : clustersPtrs){ + + int layer = HGC_layer(clu->subdetId(),clu->layer()); + if(layer>& clustersPtrs = c3d.constituents(); + std::unordered_map layers_pt; + float max_pt = 0.; + int max_layer = 0; + for(const auto& cluster_ptr : clustersPtrs){ + int layer = HGC_layer(cluster_ptr->subdetId(),cluster_ptr->layer()); + auto itr_insert = layers_pt.emplace(layer, 0.); + itr_insert.first->second += cluster_ptr->pt(); + if(itr_insert.first->second>max_pt){ + max_pt = itr_insert.first->second; + max_layer = layer; + } + } + return max_layer; +} + + +int HGCalShowerShape::lastLayer(const l1t::HGCalMulticluster& c3d) const { + + const std::vector>& clustersPtrs = c3d.constituents(); + + int lastLayer=-999; + + for(const auto& clu : clustersPtrs){ + + int layer = HGC_layer(clu->subdetId(),clu->layer()); + if(layer>lastLayer) lastLayer=layer; + + } + + return lastLayer; + +} + +int HGCalShowerShape::coreShowerLength(const l1t::HGCalMulticluster& c3d, const HGCalTriggerGeometryBase& triggerGeometry) const +{ + const std::vector>& clustersPtrs = c3d.constituents(); + std::vector layers(kLayersEE_+kLayersFH_+kLayersBH_); + for(const auto& cluster_ptr : clustersPtrs) + { + int layer = triggerGeometry.triggerLayer(cluster_ptr->detId()); + layers[layer-1] = true; + } + int length = 0; + int maxlength = 0; + for(bool layer : layers) + { + if(layer) length++; + else length = 0; + if(length>maxlength) maxlength = length; + } + return maxlength; +} + + +float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const { + + const std::vector>& clustersPtrs = c3d.constituents(); + + std::vector > tc_energy_eta ; + + for(const auto& clu : clustersPtrs){ + + const std::vector>& triggerCells = clu->constituents(); + + for(const auto& tc : triggerCells){ + + tc_energy_eta.emplace_back( std::make_pair(tc->energy(),tc->eta()) ); + + } + + } + + float SeeTot = sigmaXX(tc_energy_eta,c3d.eta()); + + return SeeTot; + +} + + + + +float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const { + + const std::vector>& clustersPtrs = c3d.constituents(); + + std::vector > tc_energy_phi ; + + for(const auto& clu : clustersPtrs){ + + const std::vector>& triggerCells = clu->constituents(); + + for(const auto& tc : triggerCells){ + + tc_energy_phi.emplace_back( std::make_pair(tc->energy(),tc->phi()) ); + + } + + } + + float SppTot = sigmaPhiPhi(tc_energy_phi,c3d.phi()); + + return SppTot; + +} + + + + + +float HGCalShowerShape::sigmaRRTot(const l1t::HGCalMulticluster& c3d) const { + + const std::vector>& clustersPtrs = c3d.constituents(); + + std::vector > tc_energy_r ; + + for(const auto& clu : clustersPtrs){ + + const std::vector>& triggerCells = clu->constituents(); + + for(const auto& tc : triggerCells){ + + float r = std::sqrt( pow(tc->position().x(),2) + pow(tc->position().y(),2) )/std::abs(tc->position().z()); + tc_energy_r.emplace_back( std::make_pair(tc->energy(),r) ); + + } + + } + + float r_mean = meanX(tc_energy_r); + float Szz = sigmaXX(tc_energy_r,r_mean); + + return Szz; + +} + + + + + +float HGCalShowerShape::sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const { + + std::unordered_map > > tc_layer_energy_eta; + std::unordered_map layer_LV; + + const std::vector>& clustersPtrs = c3d.constituents(); + + for(const auto& clu : clustersPtrs){ + + int layer = HGC_layer(clu->subdetId(),clu->layer()); + + layer_LV[layer] += clu->p4(); + + const std::vector>& triggerCells = clu->constituents(); + + for(const auto& tc : triggerCells){ + + tc_layer_energy_eta[layer].emplace_back( std::make_pair(tc->energy(),tc->eta()) ); + + } + + } + + + float SigmaEtaEtaMax=0; + + for(auto& tc_iter : tc_layer_energy_eta){ + + const std::vector >& energy_eta_layer = tc_iter.second; + const LorentzVector& LV_layer = layer_LV[tc_iter.first]; + float SigmaEtaEtaLayer = sigmaXX(energy_eta_layer,LV_layer.eta()); //RMS wrt layer eta, not wrt c3d eta + if(SigmaEtaEtaLayer > SigmaEtaEtaMax) SigmaEtaEtaMax = SigmaEtaEtaLayer; + + } + + + return SigmaEtaEtaMax; + + +} + + + + +float HGCalShowerShape::sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const { + + std::unordered_map > > tc_layer_energy_phi; + std::unordered_map layer_LV; + + const std::vector>& clustersPtrs = c3d.constituents(); + + for(const auto& clu : clustersPtrs){ + + int layer = HGC_layer(clu->subdetId(),clu->layer()); + + layer_LV[layer] += clu->p4(); + + const std::vector>& triggerCells = clu->constituents(); + + for(const auto& tc : triggerCells){ + + tc_layer_energy_phi[layer].emplace_back( std::make_pair(tc->energy(),tc->phi()) ); + + } + + } + + + float SigmaPhiPhiMax=0; + + for(auto& tc_iter : tc_layer_energy_phi){ + + const std::vector >& energy_phi_layer = tc_iter.second; + const LorentzVector& LV_layer = layer_LV[tc_iter.first]; + float SigmaPhiPhiLayer = sigmaPhiPhi(energy_phi_layer,LV_layer.phi()); //RMS wrt layer phi, not wrt c3d phi + if(SigmaPhiPhiLayer > SigmaPhiPhiMax) SigmaPhiPhiMax = SigmaPhiPhiLayer; + + } + + + return SigmaPhiPhiMax; + + +} + + + + +float HGCalShowerShape::sigmaRRMax(const l1t::HGCalMulticluster& c3d) const { + + std::unordered_map > > tc_layer_energy_r; + + const std::vector>& clustersPtrs = c3d.constituents(); + + for(const auto& clu : clustersPtrs){ + + int layer = HGC_layer(clu->subdetId(),clu->layer()); + + const std::vector>& triggerCells = clu->constituents(); + + for(const auto& tc : triggerCells){ + + float r = std::sqrt( pow(tc->position().x(),2) + pow(tc->position().y(),2) )/std::abs(tc->position().z()); + tc_layer_energy_r[layer].emplace_back( std::make_pair(tc->energy(),r) ); + + } + + } + + + float SigmaRRMax=0; + + for(auto& tc_iter : tc_layer_energy_r){ + + const std::vector >& energy_r_layer = tc_iter.second; + float r_mean_layer = meanX(energy_r_layer); + float SigmaRRLayer = sigmaXX(energy_r_layer,r_mean_layer); + if(SigmaRRLayer > SigmaRRMax) SigmaRRMax = SigmaRRLayer; + + } + + + return SigmaRRMax; + + +} + + +float HGCalShowerShape::sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius) const { + + const std::vector>& clustersPtrs = c3d.constituents(); + // group trigger cells by layer + std::unordered_map> > layers_tcs; + for(const auto& clu : clustersPtrs){ + int layer = HGC_layer(clu->subdetId(),clu->layer()); + const std::vector>& triggerCells = clu->constituents(); + for(const auto& tc : triggerCells){ + layers_tcs[layer].emplace_back(tc); + } + } + + // Select trigger cells within X cm of the max TC in the layer + std::unordered_map > > tc_layers_energy_r; + for(const auto& layer_tcs : layers_tcs){ + int layer = layer_tcs.first; + edm::Ptr max_tc = layer_tcs.second.front(); + for(const auto& tc : layer_tcs.second){ + if(tc->energy()>max_tc->energy()) max_tc = tc; + } + for(const auto& tc : layer_tcs.second){ + double dx = tc->position().x() - max_tc->position().x(); + double dy = tc->position().y() - max_tc->position().y(); + double distance_to_max = std::sqrt(dx*dx+dy*dy); + if(distance_to_maxposition().x()*tc->position().x() + tc->position().y()*tc->position().y())/std::abs(tc->position().z()); + tc_layers_energy_r[layer].emplace_back( std::make_pair(tc->energy(), r)); + } + } + } + + // Compute srr layer by layer + std::vector> layers_energy_srr2; + for(const auto& layer_energy_r : tc_layers_energy_r){ + const auto& energies_r = layer_energy_r.second; + float r_mean_layer = meanX(energies_r); + float srr = sigmaXX(energies_r,r_mean_layer); + double energy_sum = 0.; + for(const auto& energy_r : energies_r){ + energy_sum += energy_r.first; + } + layers_energy_srr2.emplace_back(std::make_pair(energy_sum, srr*srr)); + } + // Combine all layer srr + float srr2_mean = meanX(layers_energy_srr2); + return std::sqrt(srr2_mean); +} + + + +float HGCalShowerShape::eMax(const l1t::HGCalMulticluster& c3d) const { + + std::unordered_map layer_energy; + + const std::vector>& clustersPtrs = c3d.constituents(); + + for(const auto& clu : clustersPtrs){ + + int layer = HGC_layer(clu->subdetId(),clu->layer()); + layer_energy[layer] += clu->energy(); + + } + + float EMax=0; + + for(const auto& layer : layer_energy){ + + if(layer.second>EMax) EMax = layer.second; + + } + + return EMax; + +} + + + + + +float HGCalShowerShape::sigmaZZ(const l1t::HGCalMulticluster& c3d) const { + + const std::vector>& clustersPtrs = c3d.constituents(); + + std::vector > tc_energy_z ; + + for(const auto& clu : clustersPtrs){ + + const std::vector>& triggerCells = clu->constituents(); + + for(const auto& tc : triggerCells){ + + tc_energy_z.emplace_back( std::make_pair(tc->energy(),tc->position().z()) ); + + } + + } + + float z_mean = meanX(tc_energy_z); + float Szz = sigmaXX(tc_energy_z,z_mean); + + return Szz; + +} + + + + + +float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const { + + const std::vector>& cellsPtrs = c2d.constituents(); + + std::vector > tc_energy_eta ; + + for(const auto& cell : cellsPtrs){ + + tc_energy_eta.emplace_back( std::make_pair(cell->energy(),cell->eta()) ); + + } + + float See = sigmaXX(tc_energy_eta,c2d.eta()); + + return See; + +} + + + +float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const { + + const std::vector>& cellsPtrs = c2d.constituents(); + + std::vector > tc_energy_phi ; + + for(const auto& cell : cellsPtrs){ + + tc_energy_phi.emplace_back( std::make_pair(cell->energy(),cell->phi()) ); + + } + + float Spp = sigmaPhiPhi(tc_energy_phi,c2d.phi()); + + return Spp; + +} + + + + +float HGCalShowerShape::sigmaRRTot(const l1t::HGCalCluster& c2d) const { + + const std::vector>& cellsPtrs = c2d.constituents(); + + std::vector > tc_energy_r ; + + for(const auto& cell : cellsPtrs){ + + float r = std::sqrt( pow(cell->position().x(),2) + pow(cell->position().y(),2) )/std::abs(cell->position().z()); + tc_energy_r.emplace_back( std::make_pair(cell->energy(),r) ); + + } + + float r_mean = meanX(tc_energy_r); + float Srr = sigmaXX(tc_energy_r,r_mean); + + return Srr; + +} diff --git a/L1Trigger/L1THGCal/test/testHGCalL1TGeometry_cfg.py b/L1Trigger/L1THGCal/test/testHGCalL1TGeometry_cfg.py index dd3c15593ab31..e3c38ff19fd8b 100644 --- a/L1Trigger/L1THGCal/test/testHGCalL1TGeometry_cfg.py +++ b/L1Trigger/L1THGCal/test/testHGCalL1TGeometry_cfg.py @@ -102,12 +102,6 @@ process.load('L1Trigger.L1THGCal.hgcalTriggerPrimitives_cff') # Eventually modify default geometry parameters -process.hgcalTriggerGeometryESProducer.TriggerGeometry.TriggerGeometryName = cms.string('HGCalTriggerGeometryHexLayerBasedImp1') -process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_8inch_aligned_192_432_V8_0.txt") -process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TModulesMapping = cms.FileInPath("L1Trigger/L1THGCal/data/panel_mapping_60deg_6mod_0.txt") -process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellNeighborsMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_8inch_aligned_192_432_0.txt") -process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_mapping_BH_3x3_30deg_0.txt") -process.hgcalTriggerGeometryESProducer.TriggerGeometry.L1TCellNeighborsBHMapping = cms.FileInPath("L1Trigger/L1THGCal/data/triggercell_neighbor_mapping_BH_3x3_30deg_0.txt") process.hgcaltriggergeomtester = cms.EDAnalyzer( "HGCalTriggerGeomTester" diff --git a/L1Trigger/L1THGCal/test/testHGCalL1T_RelVal_cfg.py b/L1Trigger/L1THGCal/test/testHGCalL1T_RelVal_cfg.py index f03b16feaeb40..b9d1c60fe4462 100644 --- a/L1Trigger/L1THGCal/test/testHGCalL1T_RelVal_cfg.py +++ b/L1Trigger/L1THGCal/test/testHGCalL1T_RelVal_cfg.py @@ -31,7 +31,7 @@ # Input source process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring('/store/relval/CMSSW_9_3_0_pre2/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/92X_upgrade2023_realistic_v1_2023D17noPU-v1/00000/002E1FCB-8168-E711-BD97-0CC47A4C8EA8.root') ) + fileNames = cms.untracked.vstring('/store/relval/CMSSW_9_3_0/RelValSinglePiPt25Eta1p7_2p7/GEN-SIM-DIGI-RAW/93X_upgrade2023_realistic_v2_2023D17noPU-v1/00000/240935CF-1C9B-E711-9F7D-0025905A60BE.root') ) process.options = cms.untracked.PSet( @@ -54,7 +54,6 @@ from Configuration.AlCa.GlobalTag import GlobalTag process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') - # load HGCAL TPG simulation process.load('L1Trigger.L1THGCal.hgcalTriggerPrimitives_cff') process.hgcl1tpg_step = cms.Path(process.hgcalTriggerPrimitives)