SpacePoint2TrueHitConnectorModule.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/core/FrameworkExceptions.h>
#include <framework/core/Module.h>
 
#include <framework/datastore/StoreArray.h>
#include <framework/datastore/RelationArray.h>
 
#include <tracking/spacePointCreation/SpacePoint.h>
#include <pxd/dataobjects/PXDTrueHit.h>
#include <svd/dataobjects/SVDTrueHit.h>
 
#include <string>
#include <vector>
#include <array>
#include <tuple>
 
#include <unordered_map> // needed for typedef of defaultMap
 
// root output
#include <TFile.h>
#include <TTree.h>
 
namespace Belle2 {
  struct TrueHitInfo {
 
    int m_Id; 
    double m_wU; 
    double m_wV; 
    bool m_U; 
    bool m_V; 
    TrueHitInfo() : m_Id(-1), m_wU(0.), m_wV(0.), m_U(false), m_V(false) { }
 
    explicit TrueHitInfo(int Id) : m_Id(Id), m_wU(0.), m_wV(0.), m_U(false), m_V(false) { }
 
    // /** ctor with full information */
    // TrueHitInfo(int Id, double wU, double wV, bool U, bool V) :
    //   m_Id(Id), m_wU(wU), m_wV(wV), m_U(U), m_V(V) { }
 
    // /** ctor for single Cluster SpacePoints (SVD) or PXD SpacePoints. The information which Cluster is set can be
    //  * retrieved from the SpacePoint directlly for SVD in this case!
    //  */
    // TrueHitInfo(int Id, double weight) : m_Id(Id), m_wU(weight), m_wV(0.), m_U(true), m_V(false) { }
 
    bool operator<(const TrueHitInfo& b) const
    {
      return (std::tuple<int, double>(b.m_U + b.m_V, b.m_wU + b.m_wV) < std::tuple<int, double>(m_U + m_V, m_wU + m_wV));
    }
 
    short getNClusters() const { return m_U + m_V; }
 
    void setUWeight(double weight) { m_wU = weight; m_U = true; }
 
    void setVWeight(double weight) { m_wV = weight; m_V = true; }
 
    double getWeightSum() const { return m_wU + m_wV; }
  }; // end of struct TrueHitInfo
 
  std::ostream& operator<<(std::ostream& os, const TrueHitInfo& thInfo)
  {
    os << "Id " << thInfo.m_Id << ": wU = " << thInfo.m_wU << ", setU = " << thInfo.m_U << ". wV = " <<
       thInfo.m_wV << ", setV = " << thInfo.m_V;
    return os;
  }
 
  class SpacePoint2TrueHitConnectorModule : public Module {
 
  public:
 
    SpacePoint2TrueHitConnectorModule(); 
    void initialize() override; 
    void event() override; 
    void terminate() override; 
    BELLE2_DEFINE_EXCEPTION(NoClusterToSpacePoint,
                            "Found no related Cluster for a SpacePoint!"); 
    struct RootVariables {
      std::vector<double> SpacePointULocal; 
      std::vector<double> SpacePointVLocal; 
      std::vector<double> SpacePointXGlobal; 
      std::vector<double> SpacePointYGlobal; 
      std::vector<double> SpacePointZGlobal; 
      std::vector<double> TrueHitULocal; 
      std::vector<double> TrueHitVLocal; 
      std::vector<double> TrueHitXGlobal; 
      std::vector<double> TrueHitYGlobal; 
      std::vector<double> TrueHitZGlobal; 
      std::vector<double> WeightU; 
      std::vector<double> WeightV; 
      std::vector<unsigned short int> HitVxdID; 
      std::vector<unsigned short int> RelationStatus; 
      std::vector<unsigned int> NRelations; 
      // TODO TODO TODO TODO TODO TODO TODO: remove if not needed, only for tessting at the moment (i.e. do not commit)
      std::vector<unsigned int> ClusterSizeU; 
      std::vector<unsigned int> ClusterSizeV; 
      std::vector<double> SpacePointErrorU; 
      std::vector<double> SpacePointErrorV; 
      std::vector<double> SpacePointErrorX; 
      std::vector<double> SpacePointErrorY; 
      std::vector<double> SpacePointErrorZ; 
      // TODO TODO TODO TODO TODO TODO TODO: remove if not needed, only for tessting at the moment (i.e. do not commit)
    }; // end RootVariables
 
  protected:
 
    typedef std::unordered_map<int, TrueHitInfo> baseMapT; 
    enum e_detTypes {
      c_PXD = Belle2::VXD::SensorInfoBase::PXD, 
      c_SVD = Belle2::VXD::SensorInfoBase::SVD, 
    };
 
    enum e_relationStatus {
      c_clearHit = 1, //< bit 0: Only one TrueHit to SpacePoint
      c_ghostHit = 2, //< bit 1: ghost hit as described in this module
      c_SpacePointU = 4, //< bit 2: SpacePoint U-coordinate is set
      c_SpacePointV = 8, //< bit 3: SpacePoint V-coordinate is set
      c_primaryParticle = 16, //< bit 4: Particle related to TrueHit is primary
      c_bgParticle = 32, //< bit 5: particle related to TrueHit is background
      c_nonUniqueRelation = 64, //< bit 6: more than one TrueHit related to SpacePoint
      c_registeredRelation = 128, //< bit 7: this relation got registered
      c_noiseCluster = 256, //< bit 8: one of the Clusters has no relation to a TrueHit (i.e. is noise) -> only for twoCluster SPs!
    }; // end e_relationStatus
 
    // ================================================== PARAMETERS ==============================================================
    std::string m_PARAMoutputSuffix; 
    std::vector<std::string> m_PARAMtrueHitNames; 
    std::vector<std::string> m_PARAMspacePointNames; 
    std::vector<std::string> m_PARAMdetectorTypes;
 
    std::vector<std::string> m_PARAMclusterNames; 
    std::vector<std::string> m_PARAMrootFileName; 
    bool m_PARAMstoreSeparate;
 
    bool m_PARAMregisterAll;
 
    bool m_PARAMpositionAnalysis; 
    bool m_PARAMrequirePrimary; 
    bool m_PARAMrequireProximity; 
    double m_PARAMmaxGlobalDiff; 
//     double m_PARAMmaxLocalDiff; /**< maximum difference of local position coordinates for each direction between TrueHit and SpacePoint */
 
    double m_PARAMmaxPosSigma; 
    double m_PARAMminWeight; 
    // ================================================= INTERMALLY USED MEMBERS ==================================================
    unsigned int m_nContainers;
 
    std::vector<e_detTypes> m_detectorTypes; 
    Belle2::StoreArray<Belle2::SVDTrueHit> m_SVDTrueHits; 
    Belle2::StoreArray<Belle2::PXDTrueHit> m_PXDTrueHits; 
    Belle2::StoreArray<Belle2::PXDCluster> m_PXDClusters; 
    Belle2::StoreArray<Belle2::SVDCluster> m_SVDClusters; 
    std::vector<std::pair<Belle2::StoreArray<Belle2::SpacePoint>, e_detTypes> > m_inputSpacePoints;
 
    std::vector<Belle2::StoreArray<Belle2::SpacePoint> > m_outputSpacePoints; 
    double m_maxGlobalDiff; 
    unsigned int m_iCont;
 
    RootVariables m_rootVariables; 
    TFile* m_rootFilePtr; 
    TTree* m_treePtr; 
    // ================================================= COUNTERS =================================================================
    std::vector<unsigned int> m_SpacePointsCtr; 
    std::vector<std::array<unsigned int, 5> > m_nRelTrueHitsCtr;
 
    std::vector<unsigned int> m_noClusterCtr;
 
    std::vector<unsigned int> m_regRelationsCtr; 
    std::vector<unsigned int> m_ghostHitCtr; 
    std::vector<unsigned int> m_noTrueHitCtr;
 
    std::vector<unsigned int> m_rejectedRelsCtr;
 
    unsigned int m_weightTooSmallCtr; 
    unsigned int m_rejectedNoPrimaryCtr; 
    // TODO: make these counters for every container not only for all together!
//     unsigned int m_negWeightCtr; /**< number of negative weights */
 
//     unsigned int m_totWeightsCtr; /**< total number of weights */
 
//     unsigned int m_single2WTHCtr; /**< counter for SpacePoints with more than two possible TrueHits, but only one of them has two weights */
//
//     unsigned int m_nonSingle2WTHCtr; /**< counter for SpacePoints with more than two possible TrueHits, with more than one of them having two weights */
//
//     unsigned int m_all2WTHCtr; /**< Counter for SpacePoints where alle possible TrueHits have two weights */
//
//     unsigned int m_accSingle2WTHCtr; /**< counter for SpacePoints with more than two possible TrueHits, but only one of them has two weights, where a relation was registered to a TrueHit */
//
//     unsigned int m_accNonSingle2WTHCtr; /**< counter for SpacePoints with more than two possible TrueHits, with more than one of them having two weights where a relation was registered to a TrueHit */
//
//     unsigned int m_accAll2WTHCtr; /**< Counter for SpacePoints where alle possible TrueHits have two weights, where a relation was registered to a TrueHit */
//
//     unsigned int m_oneCluster2THCtr; /**< Counter for SpacePoints with only one Cluster but two possible TrueHits (in these cases, the one with the bigger weight gets automatically accepted) */
//
//     unsigned int m_moreThan2Weights; /**< Count the cases whith more than two weights */
//
 
    void initializeCounters();
 
    template <typename MapType, typename TrueHitType>
    std::pair<TrueHitType*, double> getTHwithWeight(const MapType& aMap, Belle2::StoreArray<TrueHitType> trueHits,
                                                    Belle2::SpacePoint* spacePoint, e_detTypes detType);
 
    template <typename TrueHitType>
    bool compatibleCombination(Belle2::SpacePoint* spacePoint, TrueHitType* trueHit);
 
    template<typename MapType, typename ClusterType, typename TrueHitType>
    MapType getRelatedTrueHits(Belle2::SpacePoint* spacePoint, std::string clusterName = "ALL", std::string trueHitName = "ALL");
 
    template<typename MapType>
    MapType processSpacePoint(Belle2::SpacePoint* spacePoint, e_detTypes detType);
 
    template<typename MapType>
    void registerAllRelations(Belle2::SpacePoint* spacePoint, MapType trueHitMap, e_detTypes detType);
 
    template<typename TrueHitType>
    void registerOneRelation(Belle2::SpacePoint* spacePoint, std::pair<TrueHitType*, double> trueHitwWeight, e_detTypes);
 
    template<typename ClusterType>
    void
    reRegisterClusterRelations(Belle2::SpacePoint* origSpacePoint, Belle2::SpacePoint* newSpacePoint, std::string clusterName = "ALL");
 
    template<typename ClusterType>
    std::vector<std::pair<ClusterType*, double> >
    getRelatedClusters(Belle2::SpacePoint* spacePoint, std::string clusterName = "ALL");
 
    template<typename TrueHitType>
    void registerTrueHitRelation(Belle2::SpacePoint* spacePoint, int index, double weight, Belle2::StoreArray<TrueHitType> trueHits);
 
    template<typename MapType>
    void positionAnalysis(Belle2::SpacePoint* spacePoint, const MapType& trueHitMap, const int& index, e_detTypes detType);
 
    std::pair<double, double> getLocalPos(Belle2::SpacePoint* spacePoint);
 
    template<typename TrueHitType>
    std::pair<ROOT::Math::XYZVector, ROOT::Math::XYZVector> getTrueHitPositions(TrueHitType* trueHit);
 
    void initializeRootFile();
 
    void closeRootFile();
 
    // TODO TODO TODO TODO TODO TODO TODO: remove if not needed, only for tessting at the moment (i.e. do not commit)
    std::pair<unsigned short int, unsigned short int> getClusterSizes(Belle2::SpacePoint* spacePoint, e_detTypes detType);
 
    std::pair<double, double> getLocalError(Belle2::SpacePoint* spacePoint);
 
    double calculateRelationWeight(const TrueHitInfo& trueHitInfo, Belle2::SpacePoint* spacePoint);
 
  }; // end module
 
  template<typename T>
  class simpleBitfield {
 
  public:
    simpleBitfield() : __bits()  { } 
    simpleBitfield(const simpleBitfield<T>& __otherBitfield) = delete; 
    simpleBitfield<T>& operator = (simpleBitfield<T>&) = delete; 
    const T hasStatus(T __statusBits) const { return (__bits & __statusBits) == __statusBits; }
    const T getStatus() const { return __bits; }
    void setStatus(T __statusBits) { __bits = __statusBits; }
    void addStatus(T __statusBits) { __bits |= __statusBits; }
    void removeStatus(T __statusBits) { __bits &= __statusBits; }
    void clear() { __bits = T(); }
 
  private:
    T __bits; 
  }; // end class bitfield
} // end namespace Belle2