MCParticle.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/gearbox/Const.h>
#include <framework/core/FrameworkExceptions.h>
#include <framework/datastore/RelationsObject.h>
 
#include <TClonesArray.h>
#include <Math/Vector3D.h>
#include <Math/Vector4D.h>
 
#include <cmath>
#include <climits>
#include <vector>
 
 
namespace Belle2 {
  class MCParticle: public RelationsObject {
  public:
 
    //Define exceptions
    BELLE2_DEFINE_EXCEPTION(LastChildIndexOutOfRangError, "Last child index out of range!")
    BELLE2_DEFINE_EXCEPTION(NoParticleListSetError, "No Particle list set, cannot determine related particles!")
    BELLE2_DEFINE_EXCEPTION(ParticlePDGNotKnownError, "The pdg value (%1%) of the MCParticle is not known!")
 
 
    
    enum StatusBit {
      c_PrimaryParticle   = 1 << 0,
      c_StableInGenerator = 1 << 1,
      c_LeftDetector      = 1 << 2,
      c_StoppedInDetector = 1 << 3,
      c_IsVirtual         = 1 << 4,
      c_Initial           = 1 << 5,
      c_IsISRPhoton       = 1 << 6,
      c_IsFSRPhoton       = 1 << 7,
      c_IsPHOTOSPhoton    = 1 << 8,
      c_IsRadiativePhoton = c_IsISRPhoton | c_IsFSRPhoton | c_IsPHOTOSPhoton,
    };
 
 
    MCParticle() {}
 
    MCParticle(TClonesArray* plist, const MCParticle& p):
      m_plist(plist),
      m_productionTime(p.m_productionTime), m_productionVertex_x(p.m_productionVertex_x),
      m_productionVertex_y(p.m_productionVertex_y), m_productionVertex_z(p.m_productionVertex_z),
      m_decayTime(p.m_decayTime), m_decayVertex_x(p.m_decayVertex_x),
      m_decayVertex_y(p.m_decayVertex_y), m_decayVertex_z(p.m_decayVertex_z),
      m_pdg(p.m_pdg), m_mass(p.m_mass), m_energy(p.m_energy),
      m_momentum_x(p.m_momentum_x), m_momentum_y(p.m_momentum_y), m_momentum_z(p.m_momentum_z),
      m_index(p.m_index),
      m_mother(p.m_mother), m_firstDaughter(p.m_firstDaughter), m_lastDaughter(p.m_lastDaughter),
      m_secondaryPhysicsProcess(p.m_secondaryPhysicsProcess),
      m_status(p.m_status),
      m_validVertex(p.m_validVertex),
      m_seenIn(p.m_seenIn) {}
 
    int getPDG() const { return m_pdg; }
 
    unsigned int getStatus(unsigned short int bitmask = USHRT_MAX) const { return m_status & bitmask; }
 
    bool hasStatus(unsigned short int bitmask) const { return (m_status & bitmask) == bitmask; }
 
    float getMass() const { return m_mass; }
 
    float getCharge() const;
 
    float getEnergy() const { return m_energy; }
 
    bool hasValidVertex() const { return m_validVertex; }
 
    float getProductionTime() const { return m_productionTime; }
 
    float getDecayTime() const { return m_decayTime; }
 
    float getLifetime() const { return m_decayTime - m_productionTime; }
 
    ROOT::Math::XYZVector getVertex() const { return getProductionVertex(); }
 
    ROOT::Math::XYZVector getProductionVertex() const
    {
      return ROOT::Math::XYZVector(m_productionVertex_x, m_productionVertex_y, m_productionVertex_z);
    }
 
    ROOT::Math::XYZVector getMomentum() const
    {
      return ROOT::Math::XYZVector(m_momentum_x, m_momentum_y, m_momentum_z);
    }
 
    ROOT::Math::PxPyPzEVector get4Vector() const
    {
      return ROOT::Math::PxPyPzEVector(m_momentum_x, m_momentum_y, m_momentum_z, m_energy);
    }
 
 
    ROOT::Math::XYZVector getDecayVertex() const
    {
      return ROOT::Math::XYZVector(m_decayVertex_x, m_decayVertex_y, m_decayVertex_z);
    }
 
    int getIndex() const { fixParticleList(); return m_index; }
 
    int getArrayIndex() const { fixParticleList(); return m_index - 1;  }
 
    int getFirstDaughter() const { return m_firstDaughter; }
 
 
    int getLastDaughter() const { return m_lastDaughter; }
 
    std::vector<Belle2::MCParticle*> getDaughters() const;
    //Need namespace qualifier because ROOT CINT has troubles otherwise
 
    const MCParticle* getDaughter(int i) const;
 
    int getNDaughters() const;
 
    MCParticle* getMother() const;
 
    int getSecondaryPhysicsProcess() const {return m_secondaryPhysicsProcess;}
 
    Const::DetectorSet getSeenInDetector() const { return m_seenIn; }
 
    bool hasSeenInDetector(Const::DetectorSet set) const { return m_seenIn.contains(set); }
 
    bool isVirtual() const;
 
    bool isInitial() const;
 
    bool isPrimaryParticle() const;
 
    void setPDG(int pdg) { m_pdg = pdg; }
 
    void setMassFromPDG();
 
    void setStatus(unsigned short int status) { m_status = status; }
 
    void addStatus(unsigned short int bitmask) { m_status |= bitmask; }
 
    void removeStatus(unsigned short int bitmask) { m_status &= (~bitmask); }
 
    void setMass(float mass) { m_mass = mass; }
 
    void setEnergy(float energy) { m_energy = energy; }
 
    void setValidVertex(bool valid) { m_validVertex = valid; }
 
    void setProductionTime(float time) { m_productionTime = time; }
 
    void setDecayTime(float time) { m_decayTime = time; }
 
    void setProductionVertex(const ROOT::Math::XYZVector& vertex)
    {
      m_productionVertex_x = vertex.X(); m_productionVertex_y = vertex.Y(), m_productionVertex_z = vertex.Z();
    }
 
    void setProductionVertex(float x, float y, float z)
    {
      m_productionVertex_x = x; m_productionVertex_y = y; m_productionVertex_z = z;
    }
 
    void setMomentum(const ROOT::Math::XYZVector& momentum)
    {
      m_momentum_x = momentum.X(); m_momentum_y = momentum.Y(), m_momentum_z = momentum.Z();
    }
 
    void setMomentum(float px, float py, float pz)
    {
      m_momentum_x = px, m_momentum_y = py; m_momentum_z = pz;
    }
 
    void set4Vector(const ROOT::Math::PxPyPzEVector& p4)
    {
      m_momentum_x = p4.px(); m_momentum_y = p4.py(); m_momentum_z = p4.pz(); m_energy = p4.energy();
    }
 
    void setDecayVertex(const ROOT::Math::XYZVector& vertex)
    {
      m_decayVertex_x = vertex.X(); m_decayVertex_y = vertex.Y(), m_decayVertex_z = vertex.Z();
    }
 
    void setDecayVertex(float x, float y, float z)
    {
      m_decayVertex_x = x; m_decayVertex_y = y; m_decayVertex_z = z;
    }
 
    void setSecondaryPhysicsProcess(int physicsProcess) { m_secondaryPhysicsProcess = physicsProcess; }
 
    void setSeenInDetector(Const::DetectorSet set) { m_seenIn = set; }
 
    void addSeenInDetector(Const::DetectorSet set) { m_seenIn += set; }
 
    void removeSeenInDetector(Const::DetectorSet set) { m_seenIn -= set; }
 
    void fixParticleList() const;
 
    void setVirtual() {  addStatus(c_IsVirtual); }
 
    void setInitial() {  addStatus(c_Initial); }
 
    virtual std::string getName() const override;
 
    virtual std::string getInfoHTML() const override;
 
    const MCParticle* getParticleFromGeneralizedIndexString(const std::string& generalizedIndex) const;
 
  protected:
 
    TClonesArray* m_plist = nullptr; 
 
    float m_productionTime = 0;     
    float m_productionVertex_x = 0; 
    float m_productionVertex_y = 0; 
    float m_productionVertex_z = 0; 
    float m_decayTime = 0;          
    float m_decayVertex_x = 0;      
    float m_decayVertex_y = 0;      
    float m_decayVertex_z = 0;      
    int m_pdg = 0;                  
    float m_mass = 0;               
    float m_energy = 0;             
    float m_momentum_x = 0;         
    float m_momentum_y = 0;         
    float m_momentum_z = 0;         
    int m_index = 0; 
 
    int m_mother = 0;               
    int m_firstDaughter = 0;        
    int m_lastDaughter = 0;         
    int m_secondaryPhysicsProcess = 0;  
    unsigned short int m_status = 0;      
    bool m_validVertex = false;         
    static const double c_epsilon;  
    Const::DetectorSet m_seenIn;  
    ClassDefOverride(MCParticle, 6); 
    friend class FixMergedObjectsModule;
  };
 
 
  inline bool MCParticle::isVirtual() const
  {
    bool virtuality = hasStatus(c_IsVirtual);
    if (!virtuality) {
      double E2 = m_energy * m_energy;
      double m2 = m_mass * m_mass;
      double p2 = m_momentum_x * m_momentum_x;
      p2 += m_momentum_y * m_momentum_y;
      p2 += m_momentum_z * m_momentum_z;
 
      virtuality = (fabs(E2 - (p2 + m2)) > c_epsilon * E2);
    }
    return virtuality;
  }
 
  inline bool MCParticle::isInitial() const
  {
    return hasStatus(c_Initial);
  }
 
  inline bool MCParticle::isPrimaryParticle() const
  {
    return hasStatus(c_PrimaryParticle);
  }
 
  inline MCParticle* MCParticle::getMother() const
  {
    if (m_mother == 0)
      return nullptr;
    fixParticleList();
    return static_cast<MCParticle*>(m_plist->At(m_mother - 1));
  }
 
} // end namespace Belle2