#include "TrackingUtils.hxx"
#include "IReconBase.hxx"
#include "IReconNode.hxx"
#include "IReconState.hxx"
#include "IHandle.hxx"
#include "ICOMETLog.hxx"
#include "IOADatabase.hxx"
#include "IFieldManager.hxx"
#include "HEPUnits.hxx"
#include "IIntegerDatum.hxx"
#include "IRealDatum.hxx"
#include "IGeomInfo.hxx"
#include "IGeomIdManager.hxx"
#include "IGeometryDatabase.hxx"

#include "IPIDManager.hxx"

#include <CLHEP/Matrix/Matrix.h>
#include <CLHEP/Matrix/Vector.h>


namespace TrackingUtils {

  //*****************************************************************************
  /// Whether an object is a IReconShower or a IReconPID that inherits from a IReconShower.
  bool IsShower(COMET::IHandle<COMET::IReconBase> object) {
  //*****************************************************************************
    COMET::IHandle<COMET::IReconShower> shower2 = object;

    if (shower2) {
      return true;
    }

    COMET::IHandle<COMET::IReconPID> pid = object;

    if (pid) {
      if (pid->GetConstituents() && pid->GetConstituents()->size() == 1) {
        COMET::IHandle<COMET::IReconShower> shower = *(pid->GetConstituents()->begin());

        if (shower) {
          return true;
        }
      }
    }

    // this is the only way shower objects will be recognised the second time they are converted
    if ( object->GetAlgorithmName()=="Shower")
      return true;

    return false;
  }

  //*****************************************************************************
  bool IsPID(const COMET::IReconBase& object){
  //*****************************************************************************

    COMET::IHandle<COMET::IReconState> tstate = object.GetState();
    if(!tstate){
      std::cout << "Strange !!! No state in object" << std::endl;
      return false;
    }

    bool ok = true;
    COMET::IHandle<COMET::IPIDState> pidstate = tstate;
    if(!pidstate) ok = false;  

    return ok;

  }

  //*****************************************************************************
  /// Get the charge for a IReconPID or IReconTrack.
  double GetCharge(COMET::IHandle<COMET::IReconBase> object) {
  //*****************************************************************************
    return GetCharge(object->GetState());
  }

  //*****************************************************************************
  /// Get the charge for a IPIDState or ITrackState.
  double GetCharge(COMET::IHandle<COMET::IReconState> state) {
  //*****************************************************************************
    // TODO: set default charge in a better way
    double defaultCharge = 1;

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) {
      double q, p;
      bool ok = Curvature_to_MomentumAndCharge(trackState->GetPosition().Vect(), trackState->GetDirection(), trackState->GetCurvature(), p, q);
      if (ok)
        return q;
      else
        return defaultCharge;
    }

    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState)
      return pidState->GetCharge();

    return 0;
  }


  //*****************************************************************************
  /// Get the momentum for a IReconPID or IReconTrack.
  double GetMomentum(COMET::IHandle<COMET::IReconBase> object) {
  //*****************************************************************************
    return GetMomentum(object->GetState());
  }

  //*****************************************************************************
  /// Get the momentum for a IPIDState or ITrackState.
  double GetMomentum(COMET::IHandle<COMET::IReconState> state) {
  //*****************************************************************************
    // TODO: try to set the default value in a better way, from RecPackConverter
    double defaultMomentum = 10000 * unit::MeV;

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) {

      double q, p;
      bool ok = Curvature_to_MomentumAndCharge(trackState->GetPosition().Vect(), trackState->GetDirection(), trackState->GetCurvature(), p, q);
      if (ok)
        return p;
      else
        return defaultMomentum;
    }

    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState)
      return pidState->GetMomentum();

    return 0;
  }

  //*****************************************************************************
  double GetMomentumError(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState){

      CLHEP::HepMatrix C = GetDirectionMomentumCov(state);    
      double ep = sqrt(C[3][3]);

      return ep;

    } 


    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState) return sqrt(pidState->GetMomentumVariance());

    return 0;

  }

  //*****************************************************************************
  CLHEP::HepMatrix GetDirectionMomentumCov(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState){

      double p = GetMomentum(state);
      double k = trackState->GetCurvature();
      double ux = trackState->GetDirection().X();

      CLHEP::HepMatrix  Jacobian(2,2,0);
      CLHEP::HepMatrix  Covariance(2,2,0);
      CLHEP::HepMatrix  Cov(2,2,0);

      Covariance[0][0] = trackState->GetCovarianceValue(trackState->GetDirectionIndex(),trackState->GetDirectionIndex());
      Covariance[0][1] = trackState->GetCovarianceValue(trackState->GetDirectionIndex(),trackState->GetCurvatureIndex());

      Covariance[1][0] = trackState->GetCovarianceValue(trackState->GetCurvatureIndex(),trackState->GetDirectionIndex());
      Covariance[1][1] = trackState->GetCovarianceValue(trackState->GetCurvatureIndex(),trackState->GetCurvatureIndex());

      Jacobian[0][0] = 1.;
      Jacobian[0][1] = 0.;

      Jacobian[1][0] = p/(1.-ux*ux)*fabs(ux);
      Jacobian[1][1] = -p/k;

      Cov = Jacobian*Covariance*Jacobian.T();


      CLHEP::HepMatrix C(4,4,1);

      C[0][3] = Cov[0][1]; 
      C[3][0] = Cov[1][0]; 
      C[3][3] = Cov[1][1]; 

      return C;
    } 


    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState) return CLHEP::HepMatrix(4,4,0);

    return CLHEP::HepMatrix(4,4,0);

  }








  //*****************************************************************************
  /// Convert the curvature of a IReconTrack into a momentum and charge. This function is copied
  /// from TRecPackConverters to avoid requirements conflicts.
  bool Curvature_to_MomentumAndCharge(const TVector3& pos, const TVector3& dir, double curv, double& p, double& q) {
  //*****************************************************************************
    double B = COMET::IFieldManager::GetFieldValue(pos).Mag() / unit::tesla;
    // project into the bending plane
    double factor = -(0.3 * B) / sqrt(1. - dir.X() * dir.X());

    if (fabs(curv) > 0 && factor != 0) {
      p = fabs(factor / curv);
      q = -curv / fabs(curv);
      return true;
    }

    return false;
  }

  //*****************************************************************************
  bool MomentumAndCharge_to_Curvature(const TVector3& pos, const TVector3& dir, double p, double q, double& curv){
  //*****************************************************************************
  
    double B = COMET::IFieldManager::GetFieldValue(pos).Mag()/unit::tesla;
    // project into the bending plane
    double factor = -(0.3*B)/sqrt(1.-dir.X()*dir.X());
    if (p > 0){
      curv = q*factor/fabs(p);
      return true;
    }
  
    return false;
  
  }



  //*****************************************************************************
  /// Set the matching_chi2ndf IDatum of an object. This datum specifies the chi-squared/NDOF for
  /// an object that has been combined from two separate objects.
  void SetMatchingChi2(COMET::IReconBase& object, double chi2ndf) {
  //*****************************************************************************
    COMET::IRealDatum * matching = new COMET::IRealDatum("matching_chi2ndf", "matching_chi2ndf");
    matching->SetValue(chi2ndf);
    object.AddDatum(matching);
  }

  //*****************************************************************************
  /// Retrieve the matching_chi2ndf IDatum from an object, if it is set.
  double GetMatchingChi2(const COMET::IReconBase& object) {
  //*****************************************************************************
    const std::string name = "matching_chi2ndf";

    if (object.Has<COMET::IRealDatum> (name))
      return object.Get<COMET::IRealDatum> (name)->GetValue();
    else
      return 0;
  }

  //*****************************************************************************
  /// Set the sense_ok IDatum of an object.
  void SetSenseOK(COMET::IReconBase& object, bool sense) {
  //*****************************************************************************
    COMET::IIntegerDatum* senseOK = new COMET::IIntegerDatum("sense_ok", "sense_ok");
    senseOK->SetValue((int) sense);
    object.AddDatum(senseOK);
  }

  //*****************************************************************************
  /// Get the sense_ok IDatum of an object.
  bool GetSenseOK(const COMET::IReconBase& object) {
  //*****************************************************************************
    const std::string name = "sense_ok";

    if (object.Has<COMET::IIntegerDatum> (name))
      return (object.Get<COMET::IIntegerDatum> (name)->GetValue() == 1);
    else
      return false;
  }

  //**************************************************************************
  void SetTRealDatumValue(COMET::IReconBase& object, double value, std::string name, std::string title){    
  //**************************************************************************

  COMET::IHandle<COMET::IRealDatum> datum =  object.Get<COMET::IRealDatum>(name.c_str());
  if(datum)
    datum->SetValue(value);
  else{
    // Add new datum
    COMET::IHandle<COMET::IRealDatum> new_datum(new COMET::IRealDatum(name.c_str(), title.c_str(), value));  
    object.AddDatum(new_datum, name.c_str());
  }

  }

  //**************************************************************************
  double GetTRealDatumValue(const COMET::IReconBase& object, std::string name, double defaultval){
    //**************************************************************************
    if (object.Has<COMET::IRealDatum>(name.c_str()))
      return object.Get<COMET::IRealDatum>(name.c_str())->GetValue();
    else
      return defaultval;
  }
  
  //******************************************************************************
  void CopyTRealDatum(COMET::IHandle<COMET::IReconBase> object_from, COMET::IHandle<COMET::IReconBase> object_to,  std::string name_from, std::string name_to){
  //******************************************************************************
    if (!object_from || !object_to) return;
    if(name_to == "")
      name_to = name_from;
  
    COMET::IHandle<COMET::IRealDatum> rdatum_from =  object_from->Get<COMET::IRealDatum>(name_from.c_str());
    if(!rdatum_from)
      return;
   
    COMET::IHandle<COMET::IRealDatum> rdatum_to = object_to->Get<COMET::IRealDatum>(name_to.c_str());
    if(rdatum_to){
      //substitute the vector
      rdatum_to->GetVector() = rdatum_from->GetVector();
    }
    else{
      COMET::IHandle<COMET::IRealDatum> new_rdatum_to(new COMET::IRealDatum(name_to.c_str(), rdatum_from->GetTitle(), 0.));
      //copy the vector
      new_rdatum_to->GetVector() = rdatum_from->GetVector();
      object_to->AddDatum(new_rdatum_to, name_to.c_str());
    }

  }

  //**************************************************************************
  void SetMergedNoRefitStatus(COMET::IReconBase& object, bool status){    
  //**************************************************************************

  COMET::IIntegerDatum* mergedNoRefit = new COMET::IIntegerDatum("merged_no_refit","merged_no_refit");
  mergedNoRefit->SetValue((bool)status);
  object.AddDatum(mergedNoRefit);    
  }

  //**************************************************************************
  bool GetMergedNoRefitStatus(const COMET::IReconBase& object){    
  //**************************************************************************

  const std::string name = "merged_no_refit";

  if(object.Has<COMET::IIntegerDatum>(name)) 
    return (object.Get<COMET::IIntegerDatum>(name)->GetValue() == 1);
  else 
    return false;
  }

  //**************************************************************************
  void SetCurvBackTrack(COMET::IReconBase& object){    
    //**************************************************************************

    COMET::IIntegerDatum* KalmanFitType = new COMET::IIntegerDatum("KalmanFitType", kKALFITCURVBACK);
    object.AddDatum(KalmanFitType);    
  }

  //**************************************************************************
  bool GetCurvBackTrack(const COMET::IReconBase& object){    
    //**************************************************************************

    bool curvBack = false;
    if (object.Has<COMET::IIntegerDatum>("KalmanFitType"))
      if (object.Get<COMET::IIntegerDatum>("KalmanFitType")->GetValue() == kKALFITCURVBACK)
        curvBack = true; 

    return curvBack;
  }









  //*****************************************************************************
  COMET::IReconPID::ParticleId GetCOMETPID(const std::string& pid1){
  //*****************************************************************************
  
  COMET::IReconPID::ParticleId pid2 = COMET::IReconPID::kNotSet; 

  if      (pid1 == "Not Set")      pid2=COMET::IReconPID::kNotSet; 
  else if (pid1 == "Other")        pid2=COMET::IReconPID::kOther; 
  else if (pid1 == "Shower")       pid2=COMET::IReconPID::kShower; 
  else if (pid1 == "EM")           pid2=COMET::IReconPID::kEM; 
  else if (pid1 == "Electron")     pid2=COMET::IReconPID::kElectron; 
  else if (pid1 == "Gamma")        pid2=COMET::IReconPID::kGamma; 
  else if (pid1 == "Hadronic")     pid2=COMET::IReconPID::kHadronic; 
  else if (pid1 == "Pi0")          pid2=COMET::IReconPID::kPiZero; 
  else if (pid1 == "LightTrack")   pid2=COMET::IReconPID::kLightTrack; 
  else if (pid1 == "Muon")         pid2=COMET::IReconPID::kMuon; 
  else if (pid1 == "Pion")         pid2=COMET::IReconPID::kPion; 
  else if (pid1 == "Heavy Track")  pid2=COMET::IReconPID::kHeavyTrack; 
  else if (pid1 == "Proton")       pid2=COMET::IReconPID::kProton; 
  else if (pid1 == "Kaon")         pid2=COMET::IReconPID::kKaon; 


  return pid2;

  }


  //*****************************************************************************
  COMET::IReconPID::ParticleId GetCOMETPID(const COMET::IReconBase& rbase){
  //*****************************************************************************

  COMET::IReconPID::ParticleId pid2 = COMET::IReconPID::kNotSet;   

  if (IsPID(rbase)){
    const COMET::IReconPID& pid = dynamic_cast<const COMET::IReconPID&>(rbase);
    pid2 = pid.GetParticleId();
  }

  return pid2;

  }

  //*****************************************************************************
  int GetNumberPID(const std::string& pid1){
  //*****************************************************************************
  
  int pid2=1;


  if      (pid1 == "Not Set")      pid2=1;
  else if (pid1 == "Other")        pid2=2;
  else if (pid1 == "Shower")       pid2=3;
  else if (pid1 == "EM")           pid2=4;
  else if (pid1 == "Electron")     pid2=5;
  else if (pid1 == "Gamma")        pid2=6;
  else if (pid1 == "Hadronic")     pid2=7;
  else if (pid1 == "Pi0")          pid2=8;
  else if (pid1 == "LightTrack")   pid2=9;
  else if (pid1 == "Muon")         pid2=10;
  else if (pid1 == "Pion")         pid2=11;
  else if (pid1 == "Heavy Track")  pid2=12;
  else if (pid1 == "Proton")       pid2=13;
  else if (pid1 == "Kaon")         pid2=14;

  return pid2;

  }

  //**************************************************************************
  UInt_t GetBit(COMET::IReconPID::ParticleId pid){    
  //**************************************************************************

  int pid1 = GetNumberPID(COMET::IReconPID::ConvertParticleId(pid));

  return BIT(pid1+1);

  }



  //**************************************************************************
  void SetParticleId(COMET::IReconState& tstate, COMET::IReconPID::ParticleId pid){    
  //**************************************************************************

  tstate.SetBit(GetBit(pid));    

  }


  //**************************************************************************
  COMET::IReconPID::ParticleId GetParticleId(const COMET::IReconState& tstate){    
  //**************************************************************************

  COMET::IReconPID::ParticleId pid2 = COMET::IReconPID::kNotSet;

  if      (tstate.TestBit(GetBit(COMET::IReconPID::kNotSet)))       pid2=COMET::IReconPID::kNotSet; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kOther)))        pid2=COMET::IReconPID::kOther; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kShower)))       pid2=COMET::IReconPID::kShower; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kEM)))           pid2=COMET::IReconPID::kEM; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kElectron)))     pid2=COMET::IReconPID::kElectron; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kGamma)))        pid2=COMET::IReconPID::kGamma; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kHadronic)))     pid2=COMET::IReconPID::kHadronic; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kPiZero)))       pid2=COMET::IReconPID::kPiZero; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kLightTrack)))   pid2=COMET::IReconPID::kLightTrack; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kMuon)))         pid2=COMET::IReconPID::kMuon; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kPion)))         pid2=COMET::IReconPID::kPion; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kHeavyTrack)))   pid2=COMET::IReconPID::kHeavyTrack; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kProton)))       pid2=COMET::IReconPID::kProton; 
  else if (tstate.TestBit(GetBit(COMET::IReconPID::kKaon)))         pid2=COMET::IReconPID::kKaon; 

  return pid2;

  }



  //*****************************************************************************
  COMET::IHandle<COMET::IReconNode> GetLastNode(const COMET::IReconBase& object, bool fitted, COMET::IReconBase::Status det){
  //*****************************************************************************


    COMET::IHandle<COMET::IReconState>  tstate;
    COMET::IHandle<COMET::IReconNode>  tnode;

    COMET::IReconNodeContainer::const_reverse_iterator it;
    for (it=object.GetNodes().rbegin(); it!=object.GetNodes().rend();it++){
      tnode = *it;
      if (det != COMET::IReconBase::kDetectorMask)
        if (!(tnode->GetObject()->UsesDetector(det))) continue;

      if (fitted){
        tstate = tnode->GetState();
        if (tstate){
          COMET::IReconPID::ParticleId pid2 = TrackingUtils::GetCOMETPID(object);
          TrackingUtils::SetParticleId(*tstate, pid2);
          return tnode;
        }
      }
      else
        return tnode;
    }

    return COMET::IHandle<COMET::IReconNode>();    

  }

  //*****************************************************************************
  COMET::IHandle<COMET::IReconNode> GetFirstNode(const COMET::IReconBase& object, bool fitted, COMET::IReconBase::Status det){
  //*****************************************************************************


    COMET::IHandle<COMET::IReconState>  tstate;
    COMET::IHandle<COMET::IReconNode>  tnode;

    COMET::IReconNodeContainer::const_iterator it;
    for (it=object.GetNodes().begin(); it!=object.GetNodes().end();it++){
      tnode = *it;
      if (det != COMET::IReconBase::kDetectorMask)
        if (!(tnode->GetObject()->UsesDetector(det))) continue;
      if (fitted){
        tstate = tnode->GetState();
        if (tstate){
          COMET::IReconPID::ParticleId pid2 = TrackingUtils::GetCOMETPID(object);
          TrackingUtils::SetParticleId(*tstate, pid2);
          return tnode;
        }
      }
      else 
        return tnode;
    }

    return COMET::IHandle<COMET::IReconNode>(); 
  }

  //*****************************************************************************
  COMET::IHandle<COMET::IReconNode> GetLastNode(const COMET::IReconBase& object, COMET::IReconBase::Status det){
  //*****************************************************************************

    return GetLastNode(object, false, det); 
  }

  //*****************************************************************************
  COMET::IHandle<COMET::IReconNode> GetFirstNode(const COMET::IReconBase& object, COMET::IReconBase::Status det){
  //*****************************************************************************

    return GetFirstNode(object, false, det); 
  }

  //*****************************************************************************
  COMET::IHandle<COMET::IReconNode> GetLastFittedNode(const COMET::IReconBase& object, COMET::IReconBase::Status det){
  //*****************************************************************************

    return GetLastNode(object, true, det);
  }

  //*****************************************************************************
  COMET::IHandle<COMET::IReconNode> GetFirstFittedNode(const COMET::IReconBase& object, COMET::IReconBase::Status det){
  //*****************************************************************************

    return GetFirstNode(object, true, det);
  }

  //*****************************************************************************
  COMET::IHandle<COMET::IReconState> GetLastState(const COMET::IReconBase& object, COMET::IReconBase::Status det){
  //*****************************************************************************

    COMET::IHandle<COMET::IReconNode>  tnode = GetLastFittedNode(object, det);
    if (tnode) return tnode->GetState();
    else return COMET::IHandle<COMET::IReconState>(); 

  }

  //*****************************************************************************
  COMET::IHandle<COMET::IReconState> GetFirstState(const COMET::IReconBase& object, COMET::IReconBase::Status det){
  //*****************************************************************************

    COMET::IHandle<COMET::IReconNode>  tnode = GetFirstFittedNode(object, det);
    if (tnode) return tnode->GetState();
    else return COMET::IHandle<COMET::IReconState>(); 

  }

  //*****************************************************************************
  TLorentzVector GetPosition(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetPosition(object->GetState());

  }

  //*****************************************************************************
  TLorentzVector GetPositionVariance(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetPositionVariance(object->GetState());

  }

  //*****************************************************************************
  TVector3 GetDirection(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetDirection(object->GetState());

  }

  //*****************************************************************************
  TVector3 GetDirectionVariance(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetDirectionVariance(object->GetState());

  }

  //*****************************************************************************
  double GetMomentumError(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetMomentumError(object->GetState());

  }


  //*****************************************************************************
  double GetEDeposit(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetEDeposit(object->GetState());

  }


  //*****************************************************************************
  TVector3 GetWidth(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetWidth(object->GetState());

  }

  //*****************************************************************************
  TVector3 GetCone(COMET::IHandle<COMET::IReconBase> object){
  //*****************************************************************************

    return GetCone(object->GetState());

  }

  //*****************************************************************************
  TLorentzVector GetPosition(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) return trackState->GetPosition();

    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState) return pidState->GetPosition();

    COMET::IHandle<COMET::IVertexState> vertexState = state;
    if (vertexState) return vertexState->GetPosition();

    COMET::IHandle<COMET::IShowerState> showerState = state;
    if (showerState) return showerState->GetPosition();

    COMET::IHandle<COMET::IClusterState> clusterState = state;
    if (clusterState) return clusterState->GetPosition();


    return TLorentzVector(0,0,0,0);

  }

  //*****************************************************************************
  TLorentzVector GetPositionVariance(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) return trackState->GetPositionVariance();

    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState) return pidState->GetPositionVariance();

    COMET::IHandle<COMET::IVertexState> vertexState = state;
    if (vertexState) return vertexState->GetPositionVariance();

    COMET::IHandle<COMET::IShowerState> showerState = state;
    if (showerState) return showerState->GetPositionVariance();

    COMET::IHandle<COMET::IClusterState> clusterState = state;
    if (clusterState) return clusterState->GetPositionVariance();


    return TLorentzVector(0,0,0,0);

  }



  //*****************************************************************************
  TVector3 GetDirection(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) return trackState->GetDirection();

    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState) return pidState->GetDirection();

    COMET::IHandle<COMET::IShowerState> showerState = state;
    if (showerState) return showerState->GetDirection();

    return TVector3(0,0,0);

  }

  //*****************************************************************************
  TVector3 GetDirectionVariance(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) return trackState->GetDirectionVariance();

    COMET::IHandle<COMET::IPIDState> pidState = state;
    if (pidState) return pidState->GetDirectionVariance();

    COMET::IHandle<COMET::IShowerState> showerState = state;
    if (showerState) return showerState->GetDirectionVariance();

    return TVector3(0,0,0);

  }












  //*****************************************************************************
  double GetEDeposit(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) return trackState->GetEDeposit();

    COMET::IHandle<COMET::IShowerState> showerState = state;
    if (showerState) return showerState->GetEDeposit();

    COMET::IHandle<COMET::IClusterState> clusterState = state;
    if (clusterState) return clusterState->GetEDeposit();


    return 0;

  }

  //*****************************************************************************
  TVector3 GetWidth(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::ITrackState> trackState = state;
    if (trackState) return trackState->GetWidth();

    return TVector3(0,0,0);

  }

  //*****************************************************************************
  TVector3 GetCone(COMET::IHandle<COMET::IReconState> state){
  //*****************************************************************************

    COMET::IHandle<COMET::IShowerState> showerState = state;
    if (showerState) return showerState->GetCone();

    return TVector3(0,0,0);

  }




  //*****************************************************************************
  COMET::IReconBase::Status GetDetector(COMET::IHandle<COMET::IHit> hit){
  //*****************************************************************************

    if (!hit) return 0;

    COMET::IGeometryId geomId = hit->GetGeomId();

    return GetDetector(geomId);

  }

  //*****************************************************************************
  COMET::IReconBase::Status GetDetector(const COMET::IGeometryId& geomId){
  //*****************************************************************************
    
    return 0;

  }


  //*****************************************************************************
  COMET::IReconBase::Status GetDetector(const TVector3& pos){
  //*****************************************************************************


    COMET::IGeometryId id; 
    bool ok = COMET::IOADatabase::Get().GeomId().GetGeometryId(pos.X(), pos.Y(), pos.Z(), id);    

    if (ok) return GetDetector(id);
    else return 0;

  }


}