Calo towers speed-up

DataFormats/CaloTowers/interface/CaloTower.h

@@ -44,10 +44,27 @@ class CaloTower : public reco::LeafCandidate {
 	    const LorentzVector& p4,
       const GlobalPoint& emPosition, const GlobalPoint& hadPosition);
 
+  CaloTower(CaloTowerDetId id, 
+	    float emE, float hadE, float outerE,
+	    int ecal_tp, int hcal_tp,
+	    GlobalVector p3, float iEnergy, bool massless,
+	    GlobalPoint emPosition, GlobalPoint hadPosition);
+
+  CaloTower(CaloTowerDetId id,
+            float emE, float hadE, float outerE,
+            int ecal_tp, int hcal_tp,
+            GlobalVector p3, float iEnergy, float imass,
+            GlobalPoint emPosition, GlobalPoint hadPosition);
+
+
 
    // setters
   void addConstituent( DetId id ) { constituents_.push_back( id ); }
   void addConstituents( const std::vector<DetId>& ids );
+#if defined(__GXX_EXPERIMENTAL_CXX0X__)
+  void setConstituents( std::vector<DetId>&& ids ) { constituents_=std::move(ids);}
+#endif
+
   void setEcalTime(int t) { ecalTime_ = t; };
   void setHcalTime(int t) { hcalTime_ = t; };
 

DataFormats/CaloTowers/src/CaloTower.cc

@@ -21,6 +21,7 @@ CaloTower::CaloTower(const CaloTowerDetId& id,
   emLvl1_(ecal_tp), hadLvl1_(hcal_tp) {}
 
 
+
 CaloTower::CaloTower(const CaloTowerDetId& id,
 		     double emE, double hadE, double outerE,
 		     int ecal_tp, int hcal_tp,
@@ -33,6 +34,29 @@ CaloTower::CaloTower(const CaloTowerDetId& id,
   emLvl1_(ecal_tp), hadLvl1_(hcal_tp) {}
 
 
+CaloTower::CaloTower(CaloTowerDetId id,
+		     float emE, float hadE, float outerE,
+		     int ecal_tp, int hcal_tp,
+		     GlobalVector p3, float iEnergy, bool massless,
+		     GlobalPoint emPos, GlobalPoint hadPos) : 
+  LeafCandidate(0, p3, iEnergy, massless, Point(0,0,0)),  
+  id_(id),
+  emPosition_(emPos), hadPosition_(hadPos),
+  emE_(emE), hadE_(hadE), outerE_(outerE),
+  emLvl1_(ecal_tp), hadLvl1_(hcal_tp) {}
+
+CaloTower::CaloTower(CaloTowerDetId id,
+                     float emE, float hadE, float outerE,
+                     int ecal_tp, int hcal_tp,
+                     GlobalVector p3, float iEnergy, float imass,
+                     GlobalPoint emPos, GlobalPoint hadPos) :
+  LeafCandidate(0, p3, iEnergy, imass, Point(0,0,0)),
+  id_(id),
+  emPosition_(emPos), hadPosition_(hadPos),
+  emE_(emE), hadE_(hadE), outerE_(outerE),
+  emLvl1_(ecal_tp), hadLvl1_(hcal_tp) {}
+
+
 // recalculated momentum-related quantities wrt user provided vertex Z position
 
 

DataFormats/Candidate/interface/LeafCandidate.h

@@ -13,6 +13,8 @@
 #include "DataFormats/Candidate/interface/iterator_imp_specific.h"
 
 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
+#include "DataFormats/GeometryVector/interface/GlobalVector.h"
+
 
 namespace reco {
 
@@ -33,6 +35,9 @@ namespace reco {
 
     typedef unsigned int index;
 
+    static double magd(GlobalVector v) { return std::sqrt(double(v.x())*double(v.x()) + double(v.y())*double(v.y()) + double(v.z())*double(v.z()) );}
+    static double dmass(GlobalVector v, double e) { double m2 = e*e-magd(v); return m2>0 ? std::sqrt(m2) : 0;}
+
     /// default constructor                                                               
     LeafCandidate() : 
       qx3_(0), pt_(0), eta_(0), phi_(0), mass_(0), 
@@ -65,19 +70,39 @@ namespace reco {
     LeafCandidate( Charge q, const LorentzVector & p4, const Point & vtx = Point( 0, 0, 0 ),
 		   int pdgId = 0, int status = 0, bool integerCharge = true ) :
       qx3_( q ), pt_( p4.pt() ), eta_( p4.eta() ), phi_( p4.phi() ), mass_( p4.mass() ),
-      vertex_( vtx ), pdgId_( pdgId ), status_( status ),
-      cachePolarFixed_( false ), cacheCartesianFixed_( false ) {
+      vertex_( vtx ), pdgId_( pdgId ), status_( status ), p4Cartesian_(p4),
+      cachePolarFixed_( false ), cacheCartesianFixed_( true ) {
       if ( integerCharge ) qx3_ *= 3;
     }
     /// constructor from values                                                           
     LeafCandidate( Charge q, const PolarLorentzVector & p4, const Point & vtx = Point( 0, 0, 0 ),
 		   int pdgId = 0, int status = 0, bool integerCharge = true ) :
       qx3_( q ), pt_( p4.pt() ), eta_( p4.eta() ), phi_( p4.phi() ), mass_( p4.mass() ),
-      vertex_( vtx ), pdgId_( pdgId ), status_( status ),
-      cachePolarFixed_( false ), cacheCartesianFixed_( false ){
+      vertex_( vtx ), pdgId_( pdgId ), status_( status ), p4Polar_(p4),
+      cachePolarFixed_( true ), cacheCartesianFixed_( false ){
       if ( integerCharge ) qx3_ *= 3;
     }
 
+    /// constructor from values  
+    LeafCandidate( Charge q, const GlobalVector & p3, float iEnergy, bool massless, const Point & vtx = Point( 0, 0, 0 ),
+		   int pdgId = 0, int status = 0, bool integerCharge = true ) :
+      qx3_( q ), pt_( p3.perp() ), eta_( p3.eta() ), phi_( p3.phi() ), mass_(massless ? 0. :  dmass(p3,iEnergy) ),
+      vertex_( vtx ), pdgId_( pdgId ), status_( status ),  p4Polar_(pt_,eta_,phi_,mass_),  p4Cartesian_(p3.x(),p3.y(),p3.z(), massless ? magd(p3) : iEnergy),
+      cachePolarFixed_( true ), cacheCartesianFixed_( true ) {
+      if ( integerCharge ) qx3_ *= 3;
+    }
+
+
+    /// constructor from values
+    LeafCandidate( Charge q, const GlobalVector & p3, float iEnergy, float imass, const Point & vtx = Point( 0, 0, 0 ),
+                   int pdgId = 0, int status = 0, bool integerCharge = true ) :
+      qx3_( q ), pt_( p3.perp() ), eta_( p3.eta() ), phi_( p3.phi() ), mass_(imass),
+      vertex_( vtx ), pdgId_( pdgId ), status_( status ),  p4Polar_(pt_,eta_,phi_,mass_),  p4Cartesian_(p3.x(),p3.y(),p3.z(), iEnergy),
+      cachePolarFixed_( true ), cacheCartesianFixed_( true ) {
+      if ( integerCharge ) qx3_ *= 3;
+    }
+
+
     /// destructor
     virtual ~LeafCandidate();
     /// first daughter const_iterator
@@ -133,7 +158,7 @@ namespace reco {
     /// energy                                                                            
     virtual double energy() const GCC11_FINAL { cacheCartesian(); return p4Cartesian_.E(); }
     /// transverse energy                                                                 
-    virtual double et() const GCC11_FINAL { cachePolar(); return p4Polar_.Et(); }
+    virtual double et() const GCC11_FINAL { cachePolar(); return (pt_<=0) ? 0 :  p4Polar_.Et(); }
     /// mass                                                                              
     virtual float mass() const GCC11_FINAL { return mass_; }
     /// mass squared                                                                      

DataFormats/Candidate/test/testCandidate.cc

@@ -3,6 +3,7 @@
 #include "DataFormats/Candidate/interface/Candidate.h"
 #include "DataFormats/Candidate/interface/CandidateWithRef.h"
 #include <memory>
+#include <iostream>
 
 class testCandidate : public CppUnit::TestFixture {
   CPPUNIT_TEST_SUITE(testCandidate);
@@ -59,7 +60,10 @@ namespace reco {
 }
 
 void testCandidate::checkAll() {
-  reco::Particle::LorentzVector p( 1.0, 2.0, 3.0, 4.0 );
+  reco::Particle::LorentzVector p( 1.0, 2.0, 3.0, 5.0 );
+  GlobalVector v(1.0, 2.0, 3.0);
+  reco::LeafCandidate::PolarLorentzVector pl(p);  
+
   reco::Particle::Charge q( 1 );
   int x = 123, y0 = 111, y1 = 222;
   std::auto_ptr<reco::Candidate> c( new test::DummyCandidate1( p, q, x, y0, y1 ) );
@@ -78,4 +82,41 @@ void testCandidate::checkAll() {
   reco::Candidate::const_iterator b = c->begin(), e = c->end();
   CPPUNIT_ASSERT( b == e );
   CPPUNIT_ASSERT( e - b == 0 );
+
+  // test constructors
+
+  reco::LeafCandidate c1(q,p);
+  reco::LeafCandidate c2(q,pl);
+  reco::LeafCandidate c3(q,v,5.f,c1.mass());
+
+  auto ftoi = [](float x)->int { int i; memcpy(&i,&x,4); return i;};
+  auto print = [](float a, float b)->bool { std::cout << "\nwhat? " << a <<' ' << b << std::endl; return false;};
+  auto ok = [&](float a, float b)->bool { return std::abs(ftoi(a)-ftoi(b))<10 ? true : print(a,b); };
+
+  CPPUNIT_ASSERT(ok(c1.pt(),c2.pt()));
+  CPPUNIT_ASSERT(ok(c1.eta(),c2.eta()));
+  CPPUNIT_ASSERT(ok(c1.phi(),c2.phi()));
+  CPPUNIT_ASSERT(ok(c1.mass(),c2.mass()));
+
+  CPPUNIT_ASSERT(ok(c1.y(),c2.y()));
+
+  CPPUNIT_ASSERT(ok(c1.pt(),c3.pt()));
+  CPPUNIT_ASSERT(ok(c1.eta(),c3.eta()));
+  CPPUNIT_ASSERT(ok(c1.phi(),c3.phi()));
+  CPPUNIT_ASSERT(ok(c1.mass(),c3.mass()));
+
+  CPPUNIT_ASSERT(ok(c1.y(),c3.y()));
+
+  CPPUNIT_ASSERT(ok(c1.px(),c2.px()));
+  CPPUNIT_ASSERT(ok(c1.py(),c2.py()));
+  CPPUNIT_ASSERT(ok(c1.pz(),c2.pz()));
+  CPPUNIT_ASSERT(ok(c1.energy(),c2.energy()));
+
+  CPPUNIT_ASSERT(ok(c1.px(),c3.px()));
+  CPPUNIT_ASSERT(ok(c1.py(),c3.py()));
+  CPPUNIT_ASSERT(ok(c1.pz(),c3.pz()));
+  CPPUNIT_ASSERT(ok(c1.energy(),c3.energy()));
+
+
+
 }

DataFormats/Common/interface/SortedCollection.h

@@ -33,6 +33,8 @@ unreliable if such duplicate entries are made.
 #include "DataFormats/Provenance/interface/ProductID.h"
 #include "FWCore/Utilities/interface/EDMException.h"
 
+#include "FWCore/Utilities/interface/GCC11Compatibility.h"
+
 #include <algorithm>
 #include <typeinfo>
 #include <vector>
@@ -99,6 +101,13 @@ namespace edm {
     //SortedCollection(InputIterator b, InputIterator e);
 
     void push_back(T const& t);
+#if defined(__GXX_EXPERIMENTAL_CXX0X__)
+    void push_back(T && t) { obj.push_back(t);}
+
+    template<typename... Args >
+    void emplace_back( Args&&... args ) { obj.emplace_back(args...);}
+#endif
+    void pop_back() { obj.pop_back(); }
 
     void swap(SortedCollection& other);
 

Geometry/CaloGeometry/interface/CaloCellGeometry.h

@@ -73,7 +73,9 @@ class CaloCellGeometry
   virtual const CornersVec& getCorners() const = 0 ;
 
   /// Returns the position of reference for this cell 
-  const GlobalPoint& getPosition() const {return m_refPoint ; }
+  const GlobalPoint& getPosition() const {return m_refPoint;}
+  const GlobalPoint& getBackPoint() const {return m_backPoint;} 
+
   float etaPos() const { return m_eta;}
   float phiPos() const { return m_phi;}
 
@@ -123,10 +125,21 @@ class CaloCellGeometry
                       getCorners()[2].eta());
      m_dPhi = std::abs(getCorners()[0].phi() -
                       getCorners()[2].phi());
+     initBack();
+
+  }
+
+ void initBack() const {
+    // from CaloTower code
+    CornersVec const & cv = getCorners();
+    m_backPoint = GlobalPoint(0.25 * (cv[4].x() + cv[5].x() + cv[6].x() + cv[7].x()),
+                              0.25 * (cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y()),
+                              0.25 * (cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z()));   
   }
 
 private:
   GlobalPoint         m_refPoint ;
+  mutable GlobalPoint         m_backPoint ;
   mutable CornersVec  m_corners  ;
   const CCGFloat*     m_parms    ;
   float m_eta, m_phi;

Geometry/CaloTopology/src/CaloTowerConstituentsMap.cc

@@ -60,6 +60,10 @@ void CaloTowerConstituentsMap::assign(const DetId& cell, const CaloTowerDetId& t
 
 void CaloTowerConstituentsMap::sort() {
   m_items.sort();
+
+//  for (auto const & it : m_items)
+//    std::cout << std::hex << it.cell.rawId() << " " << it.tower.rawId() << std::dec << std::endl;
+
 }
 
 std::vector<DetId> CaloTowerConstituentsMap::constituentsOf(const CaloTowerDetId& id) const {

RecoJets/JetProducers/plugins/VirtualJetProducer.cc

@@ -410,6 +410,7 @@ void VirtualJetProducer::inputTowers( )
     inEnd = inputs_.end(), i = inBegin;
   for (; i != inEnd; ++i ) {
     reco::CandidatePtr input = *i;
+    // std::cout << "CaloTowerVI jets " << input->pt() << " " << input->et() << ' '<< input->energy() << ' ' << (isAnomalousTower(input) ? " bad" : " ok") << std::endl; 
     if (edm::isNotFinite(input->pt()))           continue;
     if (input->et()    <inputEtMin_)  continue;
     if (input->energy()<inputEMin_)   continue;
@@ -694,6 +695,7 @@ void VirtualJetProducer::writeJets( edm::Event & iEvent, edm::EventSetup const&
 
 
     // std::cout << "area " << ijet << " " << jetArea << " " << Area<T>::get(jet) << std::endl;
+    // std::cout << "JetVI " << ijet << jet.pt() << " " << jet.et() << ' '<< jet.energy() << ' '<< jet.mass() << std::endl;
 
     // add to the list
     jets->push_back(jet);        

RecoLocalCalo/CaloTowersCreator/interface/CaloTowersCreationAlgo.h

@@ -23,6 +23,7 @@
 
 // need if we want to store the handles
 #include "FWCore/Framework/interface/ESHandle.h"
+#include <tuple>
 
 
 #include <map>
@@ -46,6 +47,9 @@ class DetId;
 
 class CaloTowersCreationAlgo {
 public:
+
+  int nalgo=-1;
+
   CaloTowersCreationAlgo();
 
   CaloTowersCreationAlgo(double EBthreshold, double EEthreshold, 
@@ -111,6 +115,8 @@ class CaloTowersCreationAlgo {
   // The key is the calotower id.
   void makeHcalDropChMap();
 
+  void makeEcalBadChs();
+
   void begin();
   void process(const HBHERecHitCollection& hbhe);
   void process(const HORecHitCollection& ho);
@@ -139,7 +145,8 @@ class CaloTowersCreationAlgo {
   // Called in assignHit to check if the energy should be added to
   // calotower, and how to flag the channel
   unsigned int hcalChanStatusForCaloTower(const CaloRecHit* hit);
-  unsigned int ecalChanStatusForCaloTower(const CaloRecHit* hit);
+
+  std::tuple<unsigned int,bool> ecalChanStatusForCaloTower(const CaloRecHit* hit);
 
   // Channel flagging is based on acceptable severity levels specified in the
   // configuration file. These methods are used to pass the values read in
@@ -192,31 +199,33 @@ class CaloTowersCreationAlgo {
   GlobalPoint hadSegmentShwrPos(DetId detId, float fracDepth);
   // "effective" point for the EM/HAD shower in CaloTower
   //  position based on non-zero energy cells
-  GlobalPoint hadShwrPos(const std::vector<std::pair<DetId,double> >& metaContains,
+  GlobalPoint hadShwrPos(const std::vector<std::pair<DetId,float> >& metaContains,
     float fracDepth, double hadE);
-  GlobalPoint emShwrPos(const std::vector<std::pair<DetId,double> >& metaContains, 
+  GlobalPoint emShwrPos(const std::vector<std::pair<DetId,float> >& metaContains, 
     float fracDepth, double totEmE);
 
   // overloaded function to get had position based on all had cells in the tower
   GlobalPoint hadShwrPos(CaloTowerDetId id, float fracDepth);
   GlobalPoint hadShwPosFromCells(DetId frontCell, DetId backCell, float fracDepth);
 
   // for Chris
-  GlobalPoint emShwrLogWeightPos(const std::vector<std::pair<DetId,double> >& metaContains, 
+  GlobalPoint emShwrLogWeightPos(const std::vector<std::pair<DetId,float> >& metaContains, 
     float fracDepth, double totEmE);
 
 
 private:
 
   struct MetaTower {
-    MetaTower();
-    double E, E_em, E_had, E_outer;
+    MetaTower(){}
+    bool empty() const { return metaConstituents.empty();}
     // contains also energy of RecHit
-    std::vector< std::pair<DetId, double> > metaConstituents;
-    double emSumTimeTimesE, hadSumTimeTimesE, emSumEForTime, hadSumEForTime; // Sum(Energy x Timing) : intermediate container
+    std::vector< std::pair<DetId, float> > metaConstituents;
+    CaloTowerDetId id;
+    float E=0, E_em=0, E_had=0, E_outer=0;
+    float emSumTimeTimesE=0, hadSumTimeTimesE=0, emSumEForTime=0, hadSumEForTime=0; // Sum(Energy x Timing) : intermediate container
 
     // needed to set CaloTower status word
-    int numBadEcalCells, numRecEcalCells, numProbEcalCells, numBadHcalCells, numRecHcalCells, numProbHcalCells; 
+    int numBadEcalCells=0, numRecEcalCells=0, numProbEcalCells=0, numBadHcalCells=0, numRecHcalCells=0, numProbHcalCells=0; 
 
  };
 
@@ -315,14 +324,18 @@ class CaloTowersCreationAlgo {
 
 
   // internal map
-  typedef std::map<CaloTowerDetId, MetaTower> MetaTowerMap;
+  typedef std::vector<MetaTower> MetaTowerMap;
   MetaTowerMap theTowerMap;
+  unsigned int theTowerMapSize=0;
 
   // Number of channels in the tower that were not used in RecHit production (dead/off,...).
   // These channels are added to the other "bad" channels found in the recHit collection. 
   typedef std::map<CaloTowerDetId, int> HcalDropChMap;
   HcalDropChMap hcalDropChMap;
 
+  // Number of bad Ecal channel in each tower
+  unsigned short ecalBadChs[CaloTowerDetId::kSizeForDenseIndexing];
+
   // clasification of channels in tower construction: the category definition is
   // affected by the setting in the configuration file
   //