CaloAnalysis/interface/SimCluster.h

0001 #ifndef SimDataFormats_SimCluster_h
0002 #define SimDataFormats_SimCluster_h
0003
0004 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0005 #include "DataFormats/Math/interface/LorentzVector.h"
0006 #include "DataFormats/Math/interface/Point3D.h"
0007 #include "DataFormats/Math/interface/Vector3D.h"
0008 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0009 #include "SimDataFormats/Track/interface/SimTrack.h"
0010 #include <vector>
0011
0012 #include "DataFormats/DetId/interface/DetId.h"
0013 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0014 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0015
0016 //
0017 // Forward declarations
0018 //
0019 class SimTrack;
0020 class EncodedEventId;
0021
0022 /** @brief Monte Carlo truth information used for tracking validation.
0023  *
0024  * Object with references to the original SimTrack and parent and daughter
0025  * TrackingVertices. Simulation with high (~100) pileup was taking too much
0026  * memory so the class was slimmed down and copies of the SimHits were removed.
0027  *
0028  * @author original author unknown, re-engineering and slimming by Subir Sarkar
0029  * (subir.sarkar@cern.ch), some tweaking and documentation by Mark Grimes
0030  * (mark.grimes@bristol.ac.uk).
0031  * @date original date unknown, re-engineering Jan-May 2013
0032  */
0033 class SimCluster {
0034   friend std::ostream &operator<<(std::ostream &s, SimCluster const &tp);
0035
0036 public:
0037   typedef int Charge;                                       ///< electric charge type
0038   typedef math::XYZTLorentzVectorD LorentzVector;           ///< Lorentz vector
0039   typedef math::PtEtaPhiMLorentzVector PolarLorentzVector;  ///< Lorentz vector
0040   typedef math::XYZPointD Point;                            ///< point in the space
0041   typedef math::XYZVectorD Vector;                          ///< point in the space
0042
0043   /// reference to reco::GenParticle
0044   typedef reco::GenParticleRefVector::iterator genp_iterator;
0045   typedef std::vector<SimTrack>::const_iterator g4t_iterator;
0046
0047   SimCluster();
0048
0049   SimCluster(const SimTrack &simtrk);
0050   SimCluster(EncodedEventId eventID, uint32_t particleID);  // for PU
0051
0052   // destructor
0053   ~SimCluster();
0054
0055   /** @brief PDG ID.
0056    *
0057    * Returns the PDG ID of the first associated gen particle. If there are no
0058    * gen particles associated then it returns type() from the first SimTrack. */
0059   int pdgId() const {
0060     if (genParticles_.empty())
0061       return g4Tracks_[0].type();
0062     else
0063       return (*genParticles_.begin())->pdgId();
0064   }
0065
0066   /** @brief Signal source, crossing number.
0067    *
0068    * Note this is taken from the first SimTrack only, but there shouldn't be any
0069    * SimTracks from different crossings in the SimCluster. */
0070   EncodedEventId eventId() const { return event_; }
0071
0072   uint64_t particleId() const { return particleId_; }
0073
0074   // Setters for G4 and reco::GenParticle
0075   void addGenParticle(const reco::GenParticleRef &ref) { genParticles_.push_back(ref); }
0076   void addG4Track(const SimTrack &t) { g4Tracks_.push_back(t); }
0077   /// iterators
0078   genp_iterator genParticle_begin() const { return genParticles_.begin(); }
0079   genp_iterator genParticle_end() const { return genParticles_.end(); }
0080   g4t_iterator g4Track_begin() const { return g4Tracks_.begin(); }
0081   g4t_iterator g4Track_end() const { return g4Tracks_.end(); }
0082
0083   // Getters for Embd and Sim Tracks
0084   const reco::GenParticleRefVector &genParticles() const { return genParticles_; }
0085   // Only for clusters from the signal vertex
0086   const std::vector<SimTrack> &g4Tracks() const { return g4Tracks_; }
0087
0088   /// @brief Electric charge. Note this is taken from the first SimTrack only.
0089   float charge() const { return g4Tracks_[0].charge(); }
0090   /// Gives charge in unit of quark charge (should be 3 times "charge()")
0091   int threeCharge() const { return lrintf(3.f * charge()); }
0092
0093   /// @brief Four-momentum Lorentz vector. Note this is taken from the first
0094   /// SimTrack only.
0095   const math::XYZTLorentzVectorF &p4() const { return theMomentum_; }
0096
0097   /// @brief spatial momentum vector
0098   math::XYZVectorF momentum() const { return p4().Vect(); }
0099
0100   /// @brief Vector to boost to the particle centre of mass frame.
0101   math::XYZVectorF boostToCM() const { return p4().BoostToCM(); }
0102
0103   /// @brief Magnitude of momentum vector. Note this is taken from the first
0104   /// SimTrack only.
0105   float p() const { return p4().P(); }
0106
0107   /// @brief Energy. Note this is taken from the first SimTrack only.
0108   float energy() const { return p4().E(); }
0109
0110   /// @brief Transverse energy. Note this is taken from the first SimTrack only.
0111   float et() const { return p4().Et(); }
0112
0113   /// @brief Mass. Note this is taken from the first SimTrack only.
0114   float mass() const { return p4().M(); }
0115
0116   /// @brief Mass squared. Note this is taken from the first SimTrack only.
0117   float massSqr() const { return pow(mass(), 2); }
0118
0119   /// @brief Transverse mass. Note this is taken from the first SimTrack only.
0120   float mt() const { return p4().Mt(); }
0121
0122   /// @brief Transverse mass squared. Note this is taken from the first SimTrack
0123   /// only.
0124   float mtSqr() const { return p4().Mt2(); }
0125
0126   /// @brief x coordinate of momentum vector. Note this is taken from the first
0127   /// SimTrack only.
0128   float px() const { return p4().Px(); }
0129
0130   /// @brief y coordinate of momentum vector. Note this is taken from the first
0131   /// SimTrack only.
0132   float py() const { return p4().Py(); }
0133
0134   /// @brief z coordinate of momentum vector. Note this is taken from the first
0135   /// SimTrack only.
0136   float pz() const { return p4().Pz(); }
0137
0138   /// @brief Transverse momentum. Note this is taken from the first SimTrack
0139   /// only.
0140   float pt() const { return p4().Pt(); }
0141
0142   /// @brief Momentum azimuthal angle. Note this is taken from the first
0143   /// SimTrack only.
0144   float phi() const { return p4().Phi(); }
0145
0146   /// @brief Momentum polar angle. Note this is taken from the first SimTrack
0147   /// only.
0148   float theta() const { return p4().Theta(); }
0149
0150   /// @brief Momentum pseudorapidity. Note this is taken from the simtrack
0151   /// before the calorimeter
0152   float eta() const { return p4().Eta(); }
0153
0154   /// @brief Rapidity. Note this is taken from the simtrack before the
0155   /// calorimeter
0156   float rapidity() const { return p4().Rapidity(); }
0157
0158   /// @brief Same as rapidity().
0159   float y() const { return rapidity(); }
0160
0161   /** @brief Status word.
0162    *
0163    * Returns status() from the first gen particle, or -99 if there are no gen
0164    * particles attached. */
0165   int status() const { return genParticles_.empty() ? -99 : (*genParticles_[0]).status(); }
0166
0167   static const unsigned int longLivedTag;  ///< long lived flag
0168
0169   /// is long lived?
0170   bool longLived() const { return status() & longLivedTag; }
0171
0172   /** @brief Gives the total number of SimHits, in the cluster */
0173   int numberOfSimHits() const { return nsimhits_; }
0174
0175   /** @brief Gives the total number of SimHits, in the cluster */
0176   int numberOfRecHits() const { return hits_.size(); }
0177
0178   /** @brief add rechit with fraction */
0179   void addRecHitAndFraction(uint32_t hit, float fraction) {
0180     hits_.emplace_back(hit);
0181     fractions_.emplace_back(fraction);
0182   }
0183
0184   /** @brief add rechit energy */
0185   void addHitEnergy(float energy) { energies_.emplace_back(energy); }
0186
0187   /** @brief Returns list of rechit IDs and fractions for this SimCluster */
0188   std::vector<std::pair<uint32_t, float>> hits_and_fractions() const {
0189     std::vector<std::pair<uint32_t, float>> result;
0190     for (size_t i = 0; i < hits_.size(); ++i) {
0191       result.emplace_back(hits_[i], fractions_[i]);
0192     }
0193     return result;
0194   }
0195
0196   /** @brief Returns list of rechit IDs and fractions in the barrel for this SimCluster */
0197   std::vector<std::pair<uint32_t, float>> barrel_hits_and_fractions() const {
0198     std::vector<std::pair<uint32_t, float>> result;
0199     for (size_t i = 0; i < hits_.size(); ++i) {
0200       DetId detid(hits_[i]);
0201       if (detid.subdetId() != EcalBarrel && detid.subdetId() != HcalBarrel && detid.subdetId() != HcalOuter)
0202         continue;
0203       result.emplace_back(hits_[i], fractions_[i]);
0204     }
0205     return result;
0206   }
0207
0208   /** @brief Returns list of rechit IDs and fractions in the endcap for this SimCluster */
0209   std::vector<std::pair<uint32_t, float>> endcap_hits_and_fractions() const {
0210     std::vector<std::pair<uint32_t, float>> result;
0211     for (size_t i = 0; i < hits_.size(); ++i) {
0212       DetId detid(hits_[i]);
0213       if (detid.subdetId() == EcalBarrel || detid.subdetId() == HcalBarrel || detid.subdetId() == HcalOuter)
0214         continue;
0215       result.emplace_back(hits_[i], fractions_[i]);
0216     }
0217     return result;
0218   }
0219
0220   /** @brief Returns list of rechit IDs and energies for this SimCluster */
0221   std::vector<std::pair<uint32_t, float>> hits_and_energies() const {
0222     assert(hits_.size() == energies_.size());
0223     std::vector<std::pair<uint32_t, float>> result;
0224     result.reserve(hits_.size());
0225     for (size_t i = 0; i < hits_.size(); ++i) {
0226       result.emplace_back(hits_[i], energies_[i]);
0227     }
0228     return result;
0229   }
0230
0231   /** @brief clear the hits and fractions list */
0232   void clearHitsAndFractions() {
0233     std::vector<uint32_t>().swap(hits_);
0234     std::vector<float>().swap(fractions_);
0235   }
0236
0237   /** @brief clear the energies list */
0238   void clearHitsEnergy() { std::vector<float>().swap(energies_); }
0239
0240   /** @brief returns the accumulated sim energy in the cluster */
0241   float simEnergy() const { return simhit_energy_; }
0242
0243   /** @brief add simhit's energy to cluster */
0244   void addSimHit(const PCaloHit &hit) {
0245     simhit_energy_ += hit.energy();
0246     ++nsimhits_;
0247   }
0248
0249 protected:
0250   uint64_t nsimhits_{0};
0251   EncodedEventId event_;
0252
0253   uint32_t particleId_{0};
0254   float simhit_energy_{0.f};
0255   std::vector<uint32_t> hits_;
0256   std::vector<float> fractions_;
0257   std::vector<float> energies_;
0258
0259   math::XYZTLorentzVectorF theMomentum_;
0260
0261   /// references to G4 and reco::GenParticle tracks
0262   std::vector<SimTrack> g4Tracks_;
0263   reco::GenParticleRefVector genParticles_;
0264 };
0265
0266 #endif  // SimDataFormats_SimCluster_H