TrackExtrapolateG4e.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
 
/* Basf2 headers. */
#include <framework/database/DBObjPtr.h>
#include <framework/datastore/StoreArray.h>
#include <framework/gearbox/Const.h>
#include <klm/dataobjects/bklm/BKLMElementNumbers.h>
#include <tracking/dataobjects/ExtHit.h>
 
/* Geant4 headers. */
#include <G4ErrorTrajErr.hh>
#include <G4ThreeVector.hh>
#include <G4TouchableHandle.hh>
 
/* ROOT headers. */
#include <TMatrixDSym.h>
 
/* C++ headers. */
#include <map>
#include <string>
#include <vector>
 
class G4ErrorFreeTrajState;
class G4ErrorSymMatrix;
class G4StepPoint;
class G4VPhysicalVolume;
class TVector3;
 
namespace Belle2 {
 
  class BeamPipeGeo;
  class COILGeometryPar;
  class ECLCluster;
  class EKLMElementNumbers;
  class KLMChannelStatus;
  class KLMCluster;
  class KLMElementNumbers;
  class KLMHit2d;
  class KLMLikelihoodParameters;
  class KLMMuidHit;
  class KLMMuidLikelihood;
  class KLMStripEfficiency;
  class MuidBuilder;
  class RecoTrack;
  class Track;
  class TrackClusterSeparation;
 
  namespace EKLM {
    class TransformDataGlobalAligned;
  }
 
  namespace Simulation {
    class ExtCylSurfaceTarget;
    class ExtManager;
  }

  enum VolTypes {
    VOLTYPE_CDC,
    VOLTYPE_TOP1,
    VOLTYPE_TOP2,
    VOLTYPE_TOP3,
    VOLTYPE_ARICH1,
    VOLTYPE_ARICH2,
    VOLTYPE_ARICH3,
    VOLTYPE_ECL,
    VOLTYPE_BKLM1,
    VOLTYPE_BKLM2,
    VOLTYPE_EKLM
  };

  struct ExtState {
    const Track* track;
    int pdgCode;
    bool isCosmic;
    double tof;
    double length;
    G4ThreeVector directionAtIP;
    double chi2;
    int nPoint;
    int firstBarrelLayer;
    int firstEndcapLayer;
    int lastBarrelExtLayer;
    int lastBarrelHitLayer;
    int lastEndcapExtLayer;
    int lastEndcapHitLayer;
    int extLayerPattern;
    int hitLayerPattern;
    bool escaped;
  };

  struct Intersection {
    int hit;
    bool inBarrel;
    bool isForward;
    int sector;
    int layer;
    G4ThreeVector position;
    G4ThreeVector momentum;
    G4ErrorSymMatrix covariance;
    G4ThreeVector positionAtHitPlane;
    double time;
    double chi2;
  };

  class TrackExtrapolateG4e {
 
  public:

    static TrackExtrapolateG4e* getInstance();

    ~TrackExtrapolateG4e();

    void initialize(double minPt, double minKE,
                    std::vector<Const::ChargedStable>& hypotheses);

    void initialize(double meanDt, double maxDt, double maxSeparation,
                    double maxKLMTrackClusterDistance, double maxECLTrackClusterDistance,
                    double minPt, double minKE, bool addHitsToRecoTrack, std::vector<Const::ChargedStable>& hypotheses);

    void beginRun(bool flag);

    void event(bool flag);

    void endRun(bool flag);

    void terminate(bool flag);

    void extrapolate(int pdgCode,
                     double tof,
                     const G4ThreeVector& position,
                     const G4ThreeVector& momentum,
                     const G4ErrorSymMatrix& covariance);
 
  private:

    TrackExtrapolateG4e();

    TrackExtrapolateG4e(TrackExtrapolateG4e&);

    void swim(ExtState&, G4ErrorFreeTrajState&,
              const std::vector<std::pair<ECLCluster*, G4ThreeVector> >*,
              const std::vector<std::pair<KLMCluster*, G4ThreeVector> >*,
              std::vector<std::map<const Track*, double> >*);

    void swim(ExtState&, G4ErrorFreeTrajState&);

    void findClosestTrackToKLMClusters();

    void registerVolumes();

    void getVolumeID(const G4TouchableHandle&, Const::EDetector&, int&);

    void fromG4eToPhasespace(const G4ErrorFreeTrajState&, G4ErrorSymMatrix&);

    void fromPhasespaceToG4e(const G4ThreeVector&, const G4ErrorSymMatrix&, G4ErrorTrajErr&);

    void fromPhasespaceToG4e(const TVector3&, const TMatrixDSym&, G4ErrorTrajErr&);

    ExtState getStartPoint(const Track&, int, G4ErrorFreeTrajState&);

    void createExtHit(const ExtHitStatus, const ExtState&, const G4ErrorFreeTrajState&, const G4StepPoint*, const G4TouchableHandle&);

    void createECLHit(const ExtState&, const G4ErrorFreeTrajState&, const G4StepPoint*, const G4StepPoint*, const G4TouchableHandle&,
                      const std::pair<ECLCluster*, G4ThreeVector>&, double, double);

    bool createMuidHit(ExtState&, G4ErrorFreeTrajState&, KLMMuidLikelihood*, std::vector<std::map<const Track*, double> >*);

    bool findBarrelIntersection(ExtState&, const G4ThreeVector&, Intersection&);

    bool findEndcapIntersection(ExtState&, const G4ThreeVector&, Intersection&);

    bool findMatchingBarrelHit(Intersection&, const Track*);

    bool findMatchingEndcapHit(Intersection&, const Track*);

    void adjustIntersection(Intersection&, const double*, const G4ThreeVector&, const G4ThreeVector&);

    void finishTrack(const ExtState&, KLMMuidLikelihood*, bool);

    static TrackExtrapolateG4e* m_Singleton;

    bool m_ExtInitialized;

    bool m_MuidInitialized;

    double m_MeanDt;

    double m_MaxDt;

    double m_MagneticField;

    double m_MaxDistSqInVariances;

    double m_MaxKLMTrackClusterDistance;

    double m_MaxECLTrackClusterDistance;

    double m_MinPt;

    double m_MinKE;

    Simulation::ExtManager* m_ExtMgr;

    const std::vector<Const::ChargedStable>* m_HypothesesExt;

    const std::vector<Const::ChargedStable>* m_HypothesesMuid;

    std::vector<Const::ChargedStable>* m_DefaultHypotheses;

    std::map<G4VPhysicalVolume*, enum VolTypes>* m_EnterExit;

    std::vector<G4VPhysicalVolume*>* m_BKLMVolumes;

    Simulation::ExtCylSurfaceTarget* m_TargetExt;

    Simulation::ExtCylSurfaceTarget* m_TargetMuid;

    DBObjPtr<COILGeometryPar> m_COILGeometryPar;

    DBObjPtr<BeamPipeGeo> m_BeamPipeGeo;

    double m_MinRadiusSq;

    double m_OffsetZ;

    int m_BarrelNSector;

    double m_BarrelMaxR;

    double m_BarrelMinR;

    double m_BarrelHalfLength;

    int m_OutermostActiveBarrelLayer;

    double m_BarrelPhiStripVariance[BKLMElementNumbers::getMaximalLayerNumber() + 1];

    double m_BarrelZStripVariance[BKLMElementNumbers::getMaximalLayerNumber() + 1];

    double m_BarrelScintVariance;

    double m_BarrelModuleMiddleRadius[2][BKLMElementNumbers::getMaximalSectorNumber() + 1]
    [BKLMElementNumbers::getMaximalLayerNumber() + 1];

    G4ThreeVector m_BarrelSectorPerp[BKLMElementNumbers::getMaximalSectorNumber() + 1];

    G4ThreeVector m_BarrelSectorPhi[BKLMElementNumbers::getMaximalSectorNumber() + 1];

    double m_EndcapMaxR;

    double m_EndcapMinR;

    double m_EndcapMiddleZ;

    double m_EndcapHalfLength;

    int m_OutermostActiveForwardEndcapLayer;

    int m_OutermostActiveBackwardEndcapLayer;

    double m_EndcapScintVariance;

    double m_EndcapModuleMiddleZ[BKLMElementNumbers::getMaximalLayerNumber() + 1];

    bool m_addHitsToRecoTrack = false;

    std::map<int, MuidBuilder*> m_MuidBuilderMap;

    const EKLMElementNumbers* m_eklmElementNumbers;

    const KLMElementNumbers* m_klmElementNumbers;

    const EKLM::TransformDataGlobalAligned* m_eklmTransformData;

    DBObjPtr<KLMChannelStatus> m_klmChannelStatus;

    DBObjPtr<KLMStripEfficiency> m_klmStripEfficiency;

    DBObjPtr<KLMLikelihoodParameters> m_klmLikelihoodParameters;

    StoreArray<ECLCluster> m_eclClusters;

    StoreArray<ExtHit> m_extHits;

    StoreArray<KLMHit2d> m_klmHit2ds;

    StoreArray<KLMCluster> m_klmClusters;

    StoreArray<KLMMuidHit> m_klmMuidHits;

    StoreArray<KLMMuidLikelihood> m_klmMuidLikelihoods;

    StoreArray<RecoTrack> m_recoTracks;

    StoreArray<Track> m_tracks;

    StoreArray<TrackClusterSeparation> m_trackClusterSeparations;
 
  };

} // end of namespace Belle2