ECLShower.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2010 - Belle II Collaboration                             *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Torben Ferber (ferber@physics.ubc.ca)                    *
 *               Guglielmo De Nardo (denardo@na.infn.it)                  *
 *               Alon Hershenhorn (hershen@phas.ubc.ca)                   *
 *               Poyuan Chen                                              *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#ifndef ECLSHOWER_H
#define ECLSHOWER_H
 
#include <framework/datastore/RelationsObject.h>
#include <TVector3.h>
#include <math.h>
#include <TMatrixDSym.h>
 
namespace Belle2 {
  class ECLShower : public RelationsObject {
  public:
 
    enum Hypothesis : unsigned int {
      c_muonNPhotons = 1,
      c_chargedHadron = 2,
      c_electronNPhotons = 3,
      c_nPhotons = 5,
      c_neutralHadron = 6,
      c_mergedPi0 = 7
    };
 
    enum StatusBit : unsigned int  {
      c_hasDeadCrystal = 1 << 0,
 
      c_hasHotCrystal = 1 << 1,
 
      c_hasProblematicCrystal = c_hasDeadCrystal | c_hasHotCrystal,
 
      c_hasPulseShapeDiscrimination = 1 << 2,
 
    };
 
    ECLShower()
    {
      m_isTrk = false;         
      m_status = 0;            
      m_showerId = 0;          
      m_connectedRegionId = 0; 
      m_hypothesisId = 0;      
      m_centralCellId = 0;     
      m_energy = 0.0;          
      m_energyRaw = 0.0;       
      m_theta = 0.0;           
      m_phi = 0.0;             
      m_r = 0.0;               
      m_Error[0] = 0.0;        
      m_Error[1] = 0.0;        
      m_Error[2] = 0.0;        
      m_Error[3] = 0.0;        
      m_Error[4] = 0.0;        
      m_Error[5] = 0.0;        
      m_time = 0;              
      m_deltaTime99 = 0;       
      m_energyHighestCrystal = 0.0;   
      m_lateralEnergy = 0.0;   
      m_minTrkDistance = 0.0;  
      m_trkDepth = 0.0;        
      m_showerDepth = 0.0;     
      m_numberOfCrystals = 0.0;     
      m_absZernike40 = 0.0;    
      m_absZernike51 = 0.0;    
      m_zernikeMVA = 0.0;      
      m_secondMoment = 0.0;    
      m_E1oE9 = 0.0;           
      m_E9oE21 = 0.0;          
      m_ShowerHadronIntensity = 0.0;         
      m_PulseShapeDiscriminationMVA = 0.5;        
      m_NumberOfHadronDigits = 0.0;         
      m_numberOfCrystalsForEnergy = 0.0;         
      m_nominalNumberOfCrystalsForEnergy = 0.0;         
      m_listOfCrystalsForEnergy = {}; 
    }
 
    void setIsTrack(bool val) { m_isTrk = val; }
 
    void setStatus(int Status) { m_status = Status; }
 
    void setShowerId(int ShowerId) { m_showerId = ShowerId; }
 
    void setConnectedRegionId(int connectedRegionId) { m_connectedRegionId = connectedRegionId; }
 
    void setHypothesisId(int hypothesisId) { m_hypothesisId = hypothesisId; }
 
    void setCentralCellId(int centralCellId) { m_centralCellId = centralCellId; }
 
    void setEnergy(double Energy) { m_energy = Energy; }
 
    void setEnergyRaw(double EnergySum) { m_energyRaw = EnergySum; }
 
    void setTheta(double Theta) { m_theta = Theta; }
 
    void setPhi(double Phi) { m_phi = Phi; }
 
    void setR(double R) { m_r = R; }
 
    void setCovarianceMatrix(double covArray[6])
    {
      for (unsigned int i = 0; i < 6; i++) {
        m_Error[i] = covArray[i];
      }
    }
 
    void setTime(double Time) { m_time = Time; }
 
    void setDeltaTime99(double TimeReso) { m_deltaTime99 = TimeReso; }
 
    void setEnergyHighestCrystal(double HighestEnergy) { m_energyHighestCrystal = HighestEnergy; }
 
    void setLateralEnergy(double lateralEnergy) { m_lateralEnergy = lateralEnergy; }
 
    void setMinTrkDistance(double dist) { m_minTrkDistance = dist; }
 
    void setTrkDepth(double trkDepth) { m_trkDepth = trkDepth; }
 
    void setShowerDepth(double showerDepth) { m_showerDepth = showerDepth; }
 
    void setNumberOfCrystals(double nofCrystals) { m_numberOfCrystals = nofCrystals; }
 
    void setAbsZernike40(double absZernike40) { m_absZernike40 = absZernike40; }
 
    void setAbsZernike51(double absZernike51) { m_absZernike51 = absZernike51; }
 
    void setZernikeMVA(double zernikeMVA) {m_zernikeMVA = zernikeMVA; }
 
    void setSecondMoment(double secondMoment) { m_secondMoment = secondMoment; }
 
    void setE1oE9(double E1oE9) { m_E1oE9 = E1oE9; }
 
    void setE9oE21(double E9oE21) { m_E9oE21 = E9oE21; }
 
    void setShowerHadronIntensity(double hadronIntensity) { m_ShowerHadronIntensity = hadronIntensity; }
 
    void setPulseShapeDiscriminationMVA(double mvaVal) { m_PulseShapeDiscriminationMVA = mvaVal; }
 
    void setNumberOfHadronDigits(double NumberOfHadronDigits) { m_NumberOfHadronDigits = NumberOfHadronDigits; }
 
    void setNumberOfCrystalsForEnergy(double numberOfCrystalsForEnergy) { m_numberOfCrystalsForEnergy = numberOfCrystalsForEnergy; }
 
    void setNominalNumberOfCrystalsForEnergy(double nominalNumberOfCrystalsForEnergy) { m_nominalNumberOfCrystalsForEnergy = nominalNumberOfCrystalsForEnergy; }
 
    void setListOfCrystalsForEnergy(const std::vector<unsigned int>& listofcrystals) { m_listOfCrystalsForEnergy = listofcrystals;}
 
 
    bool getIsTrack() const { return m_isTrk; }
 
    int getStatus() const { return m_status; }
 
    int getShowerId() const { return m_showerId; }
 
    int getConnectedRegionId() const { return m_connectedRegionId; }
 
    int getHypothesisId() const { return m_hypothesisId; }
 
    int getCentralCellId() const { return m_centralCellId; }
 
    double getEnergy() const { return m_energy; }
 
    double getEnergyRaw() const { return m_energyRaw; }
 
    double getTheta() const { return m_theta; }
 
    double getPhi() const { return m_phi; }
 
    double getR() const { return m_r; }
 
    void getCovarianceMatrixAsArray(double covArray[6]) const
    {
      for (unsigned int i = 0; i < 6; i++) {
        covArray[i] = m_Error[i];
      }
    }
 
    double getUncertaintyEnergy() const { return sqrt(m_Error[0]);}
 
    double getUncertaintyTheta() const { return  sqrt(m_Error[5]);}
 
    double getUncertaintyPhi() const {return sqrt(m_Error[2]);}
 
    double getTime() const { return m_time; }
 
    double getDeltaTime99() const { return m_deltaTime99; }
 
    double getEnergyHighestCrystal() const { return m_energyHighestCrystal; }
 
    double getLateralEnergy() const { return m_lateralEnergy; }
 
    double getMinTrkDistance() const { return m_minTrkDistance; }
 
    double getTrkDepth() const { return m_trkDepth; }
 
    double getShowerDepth() const { return m_showerDepth; }
 
    double getNumberOfCrystals() const { return m_numberOfCrystals; }
 
    double getAbsZernike40() const { return m_absZernike40; }
 
    double getAbsZernike51() const { return m_absZernike51; }
 
    double getZernikeMVA() const {return m_zernikeMVA; }
 
    double getSecondMoment() const { return m_secondMoment; }
 
    double getE1oE9() const { return m_E1oE9; }
 
    double getE9oE21() const { return m_E9oE21; }
 
    double getShowerHadronIntensity() const { return m_ShowerHadronIntensity; }
 
    double getPulseShapeDiscriminationMVA() const { return m_PulseShapeDiscriminationMVA; }
 
    double getNumberOfHadronDigits() const { return m_NumberOfHadronDigits; }
 
    double getNumberOfCrystalsForEnergy() const { return m_numberOfCrystalsForEnergy; }
 
    double getNominalNumberOfCrystalsForEnergy() const { return m_nominalNumberOfCrystalsForEnergy; }
 
    std::vector<unsigned int>& getListOfCrystalsForEnergy()  { return m_listOfCrystalsForEnergy; }
 
 
 
    TVector3 getMomentum() const
    {
      return TVector3(
               m_energy * sin(m_theta) * cos(m_phi),
               m_energy * sin(m_theta) * sin(m_phi),
               m_energy * cos(m_theta)
             );
    }
 
    TMatrixDSym getCovarianceMatrix3x3() const
    {
      TMatrixDSym covmatecl(3);
      covmatecl(0, 0) = m_Error[0];
      covmatecl(1, 0) = m_Error[1];
      covmatecl(1, 1) = m_Error[2];
      covmatecl(2, 0) = m_Error[3];
      covmatecl(2, 1) = m_Error[4];
      covmatecl(2, 2) = m_Error[5];
 
      //make symmetric
      for (int i = 0; i < 3; i++)
        for (int j = 0; j < i ; j++)
          covmatecl(j, i) = covmatecl(i, j);
      return covmatecl;
    }
 
    int getDetectorRegion() const
    {
      const double theta = getTheta();
 
      if (theta < 0.2164208) return 0;   // < 12.4deg
      if (theta < 0.5480334) return 1;   // < 31.4deg
      if (theta < 0.561996) return 11;   // < 32.2deg
      if (theta < 2.2462387) return 2;   // < 128.7deg
      if (theta < 2.2811453) return 13;   // < 130.7deg
      if (theta < 2.7070057) return 3;   // < 155.1deg
      else return 0;
    }
 
    bool hasStatus(unsigned short int bitmask) const { return (m_status & bitmask) == bitmask; }
 
    void addStatus(unsigned short int bitmask) { m_status |= bitmask; }
 
    bool hasHotCrystal() const;
 
    bool hasDeadCrystal() const;
 
    bool hasProblematicCrystal() const;
 
    bool hasPulseShapeDiscrimination() const;
 
    int getUniqueId() const
    {
      return 100000 * m_connectedRegionId + 1000 * m_hypothesisId + m_showerId;
    }
 
 
  private:
    bool m_isTrk;                   
    int m_status;                   
    int m_showerId;                 
    int m_connectedRegionId;        
    int m_hypothesisId;             
    int m_centralCellId;            
    Double32_t m_energy;            
    Double32_t m_energyRaw;         
    Double32_t m_theta;             
    Double32_t m_phi;               
    Double32_t m_r;                 
    Double32_t m_Error[6];          
    Double32_t m_time;              
    Double32_t m_deltaTime99;       
    Double32_t m_energyHighestCrystal;     
    Double32_t m_lateralEnergy;     
    Double32_t m_minTrkDistance;    
    Double32_t m_trkDepth;          
    Double32_t m_showerDepth;       
    Double32_t m_numberOfCrystals;       
    Double32_t m_absZernike40;      
    Double32_t m_absZernike51;      
    Double32_t m_zernikeMVA;        
    Double32_t m_secondMoment;      
    Double32_t m_E1oE9;             
    Double32_t m_E9oE21;            
    Double32_t
    m_ShowerHadronIntensity;        
    Double32_t
    m_PulseShapeDiscriminationMVA;        
    Double32_t
    m_NumberOfHadronDigits;         
    Double32_t m_numberOfCrystalsForEnergy; 
    Double32_t m_nominalNumberOfCrystalsForEnergy; 
    std::vector<unsigned int> m_listOfCrystalsForEnergy; 
    // 2: added uniqueID and highestE (TF)
    // 3: added LAT and distance to closest track and trk match flag (GDN)
    // 4: added time resolution (TF)
    // 5: clean up, float to Double32_t, and new variables (TF)
    // 6: changed names of Zernike moment variables/getters/setters to indicate that they are the absolute values of the moments (TF and AH)
    // 7: added centralCellId (TF)
    // 8: added zernikeMVA, removed absZernike20, 42, 53 (AH)
    // 9: renamed variables according to the new mdst scheme (TF)
    // 10: added getUniqueId()
    // 11: added m_ShowerHadronIntensity and m_NumberOfHadronDigits variables (SL)
    // 12: added m_PulseShapeDiscriminationMVA.  Noted m_ShowerHadronIntensity will be removed in release-04 (SL)
    // 13: made enums strongly typed
    // 14: added m_numberOfCrystalsForEnergy of crystals for energy determination
    // 15: added m_listOfCrystalsForEnergy, m_nominalNumberOfCrystalsForEnergy
    ClassDef(ECLShower, 15);
  };
 
  inline bool ECLShower::hasHotCrystal() const
  {
    return hasStatus(c_hasHotCrystal);
  }
 
  inline bool ECLShower::hasDeadCrystal() const
  {
    return hasStatus(c_hasDeadCrystal);
  }
 
  inline bool ECLShower::hasProblematicCrystal() const
  {
    return hasStatus(c_hasProblematicCrystal);
  }
 
  inline bool ECLShower::hasPulseShapeDiscrimination() const
  {
    return hasStatus(c_hasPulseShapeDiscrimination);
  }
 
} // end namespace Belle2
 
#endif