Hierarchy.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2015  Belle II Collaboration                              *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Tadeas Bilka                                             *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#pragma once
 
#include <framework/geometry/B2Vector3.h>
 
#include <genfit/StateOnPlane.h>
 
#include <root/TMatrixD.h>
#include <root/TGeoMatrix.h>
 
#include <Geant4/G4Transform3D.hh>
 
#include <map>
#include <set>
#include <vector>
 
namespace Belle2 {
  namespace alignment {
 
    typedef std::pair<unsigned short, unsigned short> DetectorLevelElement;
    typedef std::pair<std::vector<int>, TMatrixD> GlobalDerivativeSet;
    typedef std::vector<std::pair<int, double>> Constraint;
    typedef std::map<long, Constraint> Constraints;
 
    class GlobalDerivativesHierarchy {
    public:
      GlobalDerivativesHierarchy() {}
      virtual ~GlobalDerivativesHierarchy() {}
 
      void buildConstraints(Constraints& constraints);
 
      GlobalDerivativeSet buildGlobalDerivativesHierarchy(TMatrixD matrixChain, DetectorLevelElement child);
 
      template<class ChildDBObjectType, class MotherDBObjectType>
      void insert(unsigned short child, unsigned short mother, TMatrixD childToMotherParamTransform)
      {
        auto childUID = std::make_pair(ChildDBObjectType::getGlobalUniqueID(), child);
        auto parentUID = std::make_pair(MotherDBObjectType::getGlobalUniqueID(), mother);
 
        if (m_lookup.insert(std::make_pair(childUID, std::make_pair(parentUID, childToMotherParamTransform))).second) {
          // Just inserted an child element (so each child is inserted exactly once)
          // Try to insert parent with empty vect (does nothing if parent already in map)
          m_hierarchy.insert(std::make_pair(parentUID, std::vector<DetectorLevelElement>()));
          // insert child (can happen only once, so the vect is unique)
          m_hierarchy[parentUID].push_back(childUID);
          m_usedUniqueIDs.insert(ChildDBObjectType::getGlobalUniqueID());
          m_usedUniqueIDs.insert(MotherDBObjectType::getGlobalUniqueID());
        } else {
          // Update element transformation if inserted again
          m_lookup[childUID].second = childToMotherParamTransform;
        }
      }
 
      static void mergeGlobals(GlobalDerivativeSet& main, GlobalDerivativeSet additional);
 
      void printHierarchy();
 
      virtual std::vector<int> getElementLabels(DetectorLevelElement element) = 0;
 
      const std::set<unsigned short>& getUsedDBObjUniqueIDs() {return m_usedUniqueIDs;}
 
    private:
      std::pair<DetectorLevelElement, TMatrixD> getChildToMotherTransform(DetectorLevelElement child);
 
      std::map<DetectorLevelElement, std::pair<DetectorLevelElement, TMatrixD>> m_lookup;
      std::map<DetectorLevelElement, std::vector<DetectorLevelElement>> m_hierarchy;
 
      std::set<unsigned short> m_usedUniqueIDs {};
    };
 
    class LorentShiftHierarchy : public GlobalDerivativesHierarchy {
    public:
      LorentShiftHierarchy() : GlobalDerivativesHierarchy() {};
 
      std::vector<int> getElementLabels(DetectorLevelElement element) final;
 
      template<class LowestLevelDBObject>
      GlobalDerivativeSet getGlobalDerivatives(unsigned short sensor, const genfit::StateOnPlane* sop, B2Vector3D bField)
      {
        if (bField.Mag() < 1.e-10)
          return std::make_pair(std::vector<int>(), TMatrixD());
 
        DetectorLevelElement elementUID = std::make_pair(LowestLevelDBObject::getGlobalUniqueID(), sensor);
 
        GlobalDerivativeSet sensorDerivs = std::make_pair(getElementLabels(elementUID), getLorentzShiftDerivatives(sop, bField));
        mergeGlobals(sensorDerivs, buildGlobalDerivativesHierarchy(getLorentzShiftDerivatives(sop, bField), elementUID));
        return sensorDerivs;
      }
 
      TMatrixD getLorentzShiftDerivatives(const genfit::StateOnPlane* sop, B2Vector3D bField);
 
      template<class ChildDBObjectType, class MotherDBObjectType>
      void insertRelation(unsigned short child, unsigned short mother)
      {
        insert<ChildDBObjectType, MotherDBObjectType>(child, mother, TMatrixD(1, 1).UnitMatrix());
      }
 
    private:
    };
 
    class RigidBodyHierarchy : public GlobalDerivativesHierarchy {
 
    public:
 
      RigidBodyHierarchy() : GlobalDerivativesHierarchy() {}
 
      // Destructor
      ~RigidBodyHierarchy() {}
 
      std::vector<int> getElementLabels(DetectorLevelElement element) override;
 
      template<class LowestLevelDBObject>
      GlobalDerivativeSet getGlobalDerivatives(unsigned short sensor, const genfit::StateOnPlane* sop)
      {
        DetectorLevelElement elementUID = std::make_pair(LowestLevelDBObject::getGlobalUniqueID(), sensor);
 
        GlobalDerivativeSet sensorDerivs = std::make_pair(getElementLabels(elementUID), getRigidBodyDerivatives(sop));
        mergeGlobals(sensorDerivs, buildGlobalDerivativesHierarchy(getRigidBodyDerivatives(sop), elementUID));
        return sensorDerivs;
      }
 
      template<class ChildDBObjectType, class MotherDBObjectType>
      void insertG4Transform(unsigned short child, unsigned short mother, G4Transform3D childToMother)
      {
        insert<ChildDBObjectType, MotherDBObjectType>(child, mother, convertG4ToRigidBodyTransformation(childToMother));
      }
 
      template<class ChildDBObjectType, class MotherDBObjectType>
      void insertTGeoTransform(unsigned short child, unsigned short mother, TGeoHMatrix childToMother)
      {
        insert<ChildDBObjectType, MotherDBObjectType>(child, mother, convertTGeoToRigidBodyTransformation(childToMother));
      }
 
      TMatrixD getRigidBodyDerivatives(const genfit::StateOnPlane* sop);
 
      TMatrixD convertG4ToRigidBodyTransformation(G4Transform3D g4transform);
 
      TMatrixD convertTGeoToRigidBodyTransformation(TGeoHMatrix tgeo);
    private:
    };
  }
}