Particle.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <framework/datastore/RelationsObject.h>
#include <framework/gearbox/Const.h>
#include <mdst/dataobjects/ECLCluster.h>
 
#include <Math/Vector3D.h>
#include <Math/Vector4D.h>
#include <TMatrixFfwd.h>
#include <TMatrixFSymfwd.h>
 
#include <vector>
 
class TClonesArray;
 
namespace Belle2 {
  // forward declarations
  class KLMCluster;
  class Track;
  class TrackFitResult;
  class MCParticle;
  class PIDLikelihood;
  class V0;
 
  class Particle : public RelationsObject {
 
  public:
 
    enum EParticleSourceObject {
      c_Undefined     = 0,
      c_Track         = 1,
      c_ECLCluster    = 2,
      c_KLMCluster    = 3,
      c_V0            = 4,
      c_MCParticle    = 5,
      c_Composite     = 6,
      c_NoMDSTSource  = 7
    };
 
    enum EFlavorType {
      c_Unflavored = 0, 
      c_Flavored = 1, 
    };
 
    enum {c_DimPosition = 3, c_DimMomentum = 4, c_DimMatrix = 7,
          c_SizeMatrix = c_DimMatrix * (c_DimMatrix + 1) / 2
         };
 
    enum {c_Px, c_Py, c_Pz, c_E, c_X, c_Y, c_Z};
 
    enum PropertyFlags {
      c_Ordinary = 0, 
      c_IsUnspecified = 1, 
      c_IsIgnoreRadiatedPhotons = 2, 
      c_IsIgnoreIntermediate = 4, 
      c_IsIgnoreMassive = 8, 
      c_IsIgnoreNeutrino = 16, 
      c_IsIgnoreGamma = 32, 
      c_IsIgnoreBrems = 64, 
      c_IsIgnoreMisID = 128, 
      c_IsIgnoreDecayInFlight = 256, 
    };
 
    Particle();
 
    Particle(const ROOT::Math::PxPyPzEVector& momentum, const int pdgCode);
 
    Particle(const ROOT::Math::PxPyPzEVector& momentum,
             const int pdgCode,
             EFlavorType flavorType,
             const EParticleSourceObject particleType,
             const unsigned mdstIndex);
 
    Particle(const ROOT::Math::PxPyPzEVector& momentum,
             const int pdgCode,
             EFlavorType flavorType,
             const std::vector<int>& daughterIndices,
             TClonesArray* arrayPointer = nullptr);
 
    Particle(const ROOT::Math::PxPyPzEVector& momentum,
             const int pdgCode,
             EFlavorType flavorType,
             const std::vector<int>& daughterIndices,
             int properties,
             TClonesArray* arrayPointer = nullptr);
 
    Particle(const ROOT::Math::PxPyPzEVector& momentum,
             const int pdgCode,
             EFlavorType flavorType,
             const std::vector<int>& daughterIndices,
             int properties,
             const std::vector<int>& daughterProperties,
             TClonesArray* arrayPointer = nullptr);
 
    Particle(const Track* track,
             const Const::ChargedStable& chargedStable);
 
    Particle(int trackArrayIndex, const TrackFitResult* trackFit,
             const Const::ChargedStable& chargedStable);
 
    explicit Particle(const ECLCluster* eclCluster,
                      const Const::ParticleType& type = Const::photon);
 
    explicit Particle(const KLMCluster* klmCluster, const int pdgCode = Const::Klong.getPDGCode());
 
    explicit Particle(const MCParticle* MCparticle);
 
    ~Particle();
 
  public:
 
    // setters
 
    void setPDGCode(const int pdg)
    {
      m_pdgCode = pdg;
    }
 
    void set4Vector(const ROOT::Math::PxPyPzEVector& p4)
    {
      m_px = p4.Px();
      m_py = p4.Py();
      m_pz = p4.Pz();
      m_mass = p4.M();
    }
 
    void set4VectorDividingByMomentumScaling(const ROOT::Math::PxPyPzEVector& p4)
    {
      m_px = p4.Px() / m_momentumScale;
      m_py = p4.Py() / m_momentumScale;
      m_pz = p4.Pz() / m_momentumScale;
      m_mass = p4.M();
    }
 
    void setVertex(const ROOT::Math::XYZVector& vertex)
    {
      m_x = vertex.X();
      m_y = vertex.Y();
      m_z = vertex.Z();
    };
 
    void setEnergyLossCorrection(double energyLossCorrection)
    {
      m_energyLossCorrection = energyLossCorrection;
    }
 
    double getMomentumLossCorrectionFactor() const
    {
      if (m_energyLossCorrection == 0.0) {
        return 1.0;
      }
      double origP = m_momentumScale * sqrt(m_px * m_px + m_py * m_py + m_pz * m_pz);
      if (origP == 0.0) {
        return 1.0;
      }
      double origE = sqrt(origP * origP + m_mass * m_mass);
 
      double newP = sqrt((origE - m_energyLossCorrection) * (origE - m_energyLossCorrection) - m_mass * m_mass);
      return newP / origP;
    }
 
    void setMomentumScalingFactor(double momentumScalingFactor)
    {
      m_momentumScalingFactor = momentumScalingFactor;
      m_momentumScale = m_momentumScalingFactor * m_momentumSmearingFactor;
    }
 
    void setMomentumSmearingFactor(double momentumSmearingFactor)
    {
      m_momentumSmearingFactor = momentumSmearingFactor;
      m_momentumScale = m_momentumScalingFactor * m_momentumSmearingFactor;
    }
 
    void setJacobiMatrix(const TMatrixF& jacobiMatrix);
 
    void setMomentumVertexErrorMatrix(const TMatrixFSym& errMatrix);
 
    void setPValue(double pValue)
    {
      m_pValue = pValue;
    }
 
    void setProperty(const int properties)
    {
      m_properties = properties;
    }
 
    void updateMomentum(const ROOT::Math::PxPyPzEVector& p4,
                        const ROOT::Math::XYZVector& vertex,
                        const TMatrixFSym& errMatrix,
                        double pValue)
    {
      set4Vector(p4);
      setVertex(vertex);
      setMomentumVertexErrorMatrix(errMatrix);
      m_pValue = pValue;
    }
 
    void updateMass(const int pdgCode);
 
    void appendDaughter(const Particle* daughter, const bool updateType = true, const int daughterProperty = c_Ordinary);
 
    void appendDaughter(int particleIndex, const bool updateType = true)
    {
      if (updateType) {
        // is it a composite particle or fsr corrected?
        m_particleSource = c_Composite;
      }
      m_daughterIndices.push_back(particleIndex);
      m_daughterProperties.push_back(c_Ordinary);
    }
 
    void removeDaughter(const Particle* daughter, const bool updateType = true);
 
    bool replaceDaughter(const Particle* oldDaughter, Particle* newDaughter);
 
    bool replaceDaughterRecursively(const Particle* oldDaughter, Particle* newDaughter);
 
    // getters
 
    int getPDGCode(void) const
    {
      return m_pdgCode;
    }
 
    double getCharge(void) const;
 
    EFlavorType getFlavorType() const
    {
      return m_flavorType;
    }
 
    EParticleSourceObject getParticleSource() const
    {
      return m_particleSource;
    }
 
    unsigned getMdstArrayIndex(void) const
    {
      return m_mdstIndex;
    }
 
    int getProperty() const
    {
      return m_properties;
    }
 
    double getMass() const
    {
      return m_mass;
    }
 
    //double getMassError() const;
 
    double getPDGMass(void) const;
 
    double getPDGLifetime() const;
 
    double getEnergy() const
    {
      return sqrt(m_momentumScale * m_momentumScale * m_px * m_px +
                  m_momentumScale * m_momentumScale * m_py * m_py +
                  m_momentumScale * m_momentumScale * m_pz * m_pz +
                  m_mass * m_mass) - m_energyLossCorrection;
    }
 
    ROOT::Math::PxPyPzEVector get4Vector() const
    {
      double correction = getMomentumLossCorrectionFactor();
      return ROOT::Math::PxPyPzEVector(m_momentumScale * correction * m_px,
                                       m_momentumScale * correction * m_py,
                                       m_momentumScale * correction * m_pz,
                                       getEnergy());
    }
 
    ROOT::Math::XYZVector getMomentum() const
    {
      return m_momentumScale * getMomentumLossCorrectionFactor() * ROOT::Math::XYZVector(m_px, m_py, m_pz);
    };
 
    double getMomentumMagnitude() const
    {
      return getP();
    };
 
    double getP() const
    {
      return m_momentumScale * getMomentumLossCorrectionFactor() * sqrt(m_px * m_px + m_py * m_py + m_pz * m_pz);
    };
 
    double getPx() const
    {
      return m_momentumScale * getMomentumLossCorrectionFactor() * m_px;
    }
 
    double getPy() const
    {
      return m_momentumScale * getMomentumLossCorrectionFactor() * m_py;
    }
 
    double getPz() const
    {
      return m_momentumScale * getMomentumLossCorrectionFactor() * m_pz;
    }
 
    double getEffectiveMomentumScale() const
    {
      return m_momentumScale;
    }
 
    double getEnergyLossCorrection() const
    {
      return m_energyLossCorrection;
    }
    ROOT::Math::XYZVector getVertex() const
    {
      return ROOT::Math::XYZVector(m_x, m_y, m_z);
    };
 
    double getX() const
    {
      return m_x;
    }
 
    double getY() const
    {
      return m_y;
    }
 
    double getZ() const
    {
      return m_z;
    }
 
    double getPValue() const
    {
      return m_pValue;
    }
 
    TMatrixFSym getMomentumVertexErrorMatrix() const;
 
    TMatrixFSym getMomentumErrorMatrix() const;
 
    TMatrixFSym getVertexErrorMatrix() const;
 
    double getCosHelicity(const Particle* mother = nullptr) const;
 
    double getCosHelicityDaughter(unsigned iDaughter, unsigned iGrandDaughter = 0) const;
 
    double getAcoplanarity() const;
 
 
    int getMdstSource() const;
 
    unsigned getNDaughters(void) const
    {
      return m_daughterIndices.size();
    }
 
    const std::vector<int>& getDaughterIndices() const
    {
      return m_daughterIndices;
    }
 
    const std::vector<int>& getDaughterProperties() const
    {
      return m_daughterProperties;
    }
 
    const Particle* getDaughter(unsigned i) const;
 
    bool forEachDaughter(const std::function<bool(const Particle*)>& function,
                         bool recursive = true, bool includeSelf = true) const;
 
    std::vector<Belle2::Particle*> getDaughters() const;
    //Need namespace qualifier because ROOT CINT has troubles otherwise
 
    std::vector<const Belle2::Particle*> getFinalStateDaughters() const;
    //Need namespace qualifier because ROOT CINT has troubles otherwise
 
    std::vector<const Belle2::Particle*> getAllDaughters() const;
    //Need namespace qualifier because ROOT CINT has troubles otherwise
 
    std::vector<int> getMdstArrayIndices(EParticleSourceObject type) const;
 
    bool overlapsWith(const Particle* oParticle) const;
 
    bool isCopyOf(const Particle* oParticle, bool doDetailedComparison = false) const;
 
    const Track* getTrack() const;
 
    const TrackFitResult* getTrackFitResult() const;
 
    const V0* getV0() const;
 
    const PIDLikelihood* getPIDLikelihood() const;
 
    const ECLCluster* getECLCluster() const;
 
    double getECLClusterEnergy() const;
 
    const KLMCluster* getKLMCluster() const;
 
    const MCParticle* getMCParticle() const;
 
    std::string getName() const override;
 
    std::string getInfoHTML() const override;
 
    void print() const;
 
    std::vector<std::string> getExtraInfoNames() const;
 
    void removeExtraInfo();
 
    double getExtraInfo(const std::string& name) const;
 
    bool hasExtraInfo(const std::string& name) const;
 
    int getExtraInfoMap() const
    {
      return !m_extraInfo.empty() ? m_extraInfo[0] : -1;
    }
 
    unsigned int getExtraInfoSize() const
    {
      return m_extraInfo.size();
    }
 
    void writeExtraInfo(const std::string& name, const double value);
 
    void setExtraInfo(const std::string& name, double value);
 
    void addExtraInfo(const std::string& name, double value);
 
    TClonesArray* getArrayPointer() const
    {
      if (!m_arrayPointer)
        m_arrayPointer = RelationsObject::getArrayPointer();
      return m_arrayPointer;
    }
 
    int getPDGCodeUsedForFit() const
    {
      return std::abs(m_pdgCodeUsedForFit);
    }
 
    bool wasExactFitHypothesisUsed() const
    {
      return std::abs(m_pdgCodeUsedForFit) == std::abs(m_pdgCode);
    }
 
    bool isMostLikely() const;
 
    std::pair<Const::ChargedStable, const TrackFitResult*> getMostLikelyTrackFitResult() const;
 
    bool isMostLikelyTrackFitResult() const;
 
    ECLCluster::EHypothesisBit getECLClusterEHypothesisBit() const
    {
      const int pdg = abs(getPDGCode());
      if ((pdg == Const::photon.getPDGCode())
          or (pdg == Const::electron.getPDGCode())
          or (pdg == Const::muon.getPDGCode())
          or (pdg == Const::pion.getPDGCode())
          or (pdg == Const::kaon.getPDGCode())
          or (pdg == Const::proton.getPDGCode())
          or (pdg == Const::deuteron.getPDGCode())) {
        return ECLCluster::EHypothesisBit::c_nPhotons;
      } else if ((pdg == Const::Klong.getPDGCode())
                 or (pdg == Const::neutron.getPDGCode())) {
        return ECLCluster::EHypothesisBit::c_neutralHadron;
      } else {
        return ECLCluster::EHypothesisBit::c_none;
      }
    }
 
    const Particle* getParticleFromGeneralizedIndexString(const std::string& generalizedIndex) const;
 
    void updateJacobiMatrix();
 
    void fillFSPDaughters(std::vector<const Belle2::Particle*>& fspDaughters) const;
 
    void fillAllDaughters(std::vector<const Belle2::Particle*>& allDaughters) const;
 
 
  private:
 
    // persistent data members
    int m_pdgCode;  
    int m_pdgCodeUsedForFit = 0; 
    double m_mass;   
    double m_px;     
    double m_py;     
    double m_pz;     
    double m_momentumScale = 1.0; 
    double m_momentumScalingFactor = 1.0; 
    double m_momentumSmearingFactor = 1.0; 
    double m_energyLossCorrection = 0.0; 
    double m_x;      
    double m_y;      
    double m_z;      
    double m_errMatrix[c_SizeMatrix] = {}; 
    double m_jacobiMatrix[c_SizeMatrix] = {}; 
    double m_pValue;   
    std::vector<int> m_daughterIndices;  
    EFlavorType m_flavorType;  
    EParticleSourceObject m_particleSource;  
    unsigned m_mdstIndex;  
    int m_properties; 
    std::vector<int> m_daughterProperties; 
    int m_identifier = -1;
 
    std::vector<double> m_extraInfo;
 
    // transient data members
    mutable TClonesArray* m_arrayPointer; 
 
    // private methods
    void setMomentumPositionErrorMatrix(const TrackFitResult* trackFit);
 
    void resetErrorMatrix();
 
    void resetJacobiMatrix();
 
    void storeErrorMatrix(const TMatrixFSym& errMatrix);
 
    void storeJacobiMatrix(const TMatrixF& jacobiMatrix);
 
    // TODO: this can be optimized for speed
    void fillDecayChain(std::vector<int>& decayChain) const;
 
    void setFlavorType();
 
    void setMdstArrayIndex(const int arrayIndex);
 
    int generatePDGCodeFromCharge(const int chargedSign, const Const::ChargedStable& chargedStable);
 
    ClassDefOverride(Particle, 17); 
    // v8: added identifier, changed getMdstSource
    // v9: added m_pdgCodeUsedForFit
    // v10: added m_properties
    // v11: added m_daughterProperties
    // v12: renamed EParticleType m_particleType to EParticleSourceObject m_particleSource
    // v13: added m_momentumScale
    // v14: added m_jacobiMatrix
    // v15: added m_momentumScalingFactor and m_momentumSmearingFactor
    // v16: use double precision for private members
    // v17: added m_energyLossCorrection
    friend class ParticleSubset;
  };
 
} // end namespace Belle2