#include "IReconTrackerModule.hxx"

#include <functional>
#include <algorithm>

#include <IDatum.hxx>
#include <IRealDatum.hxx>
#include <IReconTrack.hxx>
#include <IReconShower.hxx>
#include <IReconNode.hxx>

#include <IGeomInfo.hxx>
#include <IG4HitSegment.hxx>
#include <IMultiHit.hxx>
#include <IComboHit.hxx>
#include <retrieveHitTruthInfo.hxx>


#include <IOARuntimeParameters.hxx>


#include "IRecPackManager.hxx"
#include "IValidationUtils.hxx"
#include "IEventFolder.hxx"

#include <IFieldManager.hxx>
#include <TrackTruthInfo.hxx>

ClassImp(COMET::IReconTrackerModule);
ClassImp(COMET::IReconTrackerModule::IUnusedHit);
COMET::IReconTrackerModule::IUnusedHit::~IUnusedHit(){}

ClassImp(COMET::IReconTrackerModule::ITrackerVertex);
COMET::IReconTrackerModule::ITrackerVertex::~ITrackerVertex(){}

ClassImp(COMET::IReconTrackerModule::ITPCTrack);
COMET::IReconTrackerModule::ITPCTrack::~ITPCTrack(){}

ClassImp(COMET::IReconTrackerModule::IFGDTrack);
COMET::IReconTrackerModule::IFGDTrack::~IFGDTrack(){}

ClassImp(COMET::IReconTrackerModule::ITrueParticle);
COMET::IReconTrackerModule::ITrueParticle::~ITrueParticle(){}

ClassImp(COMET::IReconTrackerModule::ITrackerNode);
COMET::IReconTrackerModule::ITrackerNode::~ITrackerNode(){}

ClassImp(COMET::IReconTrackerModule::ITrackerConstituent);
COMET::IReconTrackerModule::ITrackerConstituent::~ITrackerConstituent(){}
COMET::IReconTrackerModule::ITrackerConstituent::ITrackerConstituent(){
  NConstituents = 0;
  return;
}

ClassImp(COMET::IReconTrackerModule::ITrackerResult);
COMET::IReconTrackerModule::ITrackerResult::~ITrackerResult(){}
COMET::IReconTrackerModule::ITrackerResult::ITrackerResult(){
  Constituents = new TClonesArray("COMET::IReconTrackerModule::ITrackerConstituent", 10);
  TPC = new TClonesArray("COMET::IReconTrackerModule::ITPCTrack", 10);
  FGD = new TClonesArray("COMET::IReconTrackerModule::IFGDTrack", 10);
  Nodes = new TClonesArray("COMET::IReconTrackerModule::ITrackerNode",10);

  NConstituents = 0;
  NTotalConstituents = 0;
  Constituents->Clear();

  NTPCs = 0;
  TPC->Clear();
  
  NFGDs = 0;
  FGD->Clear();

  NNodes = 0;
  Nodes->Clear();
  
  return;
}

ClassImp(COMET::IReconTrackerModule::ITrackOther);
COMET::IReconTrackerModule::ITrackOther::~ITrackOther(){}
COMET::IReconTrackerModule::ITrackOther::ITrackOther(){
  NHits = 0;
  Hits = new TClonesArray("COMET::IReconTrackerModule::IUnusedHit", 200);
  return;
}


#define CVSTAG "\
  $Name: v5r19 $"
#define CVSID  "\
"  

//*************************************************************
COMET::IReconTrackerModule::IReconTrackerModule(const char *name, const char *title)
{
//*************************************************************
  SetNameTitle(name, title);
  // Enable this module by default:
  fIsEnabled = kTRUE;
  fDescription = "Tracker Recon Output";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;
  // Tree variables
  fNVertices = 0;
  fVertices = new TClonesArray("COMET::IReconTrackerModule::ITrackerVertex", 200);
  fNTracks = 0;
  fTracks = new TClonesArray("COMET::IReconTrackerModule::ITrackerResult", 200);
  fNFGDOther =0;
  fFGDOther = new TClonesArray("COMET::IReconTrackerModule::ITrackOther", 200);
  fNTPCOther =0;
  fTPCOther = new TClonesArray("COMET::IReconTrackerModule::ITrackOther", 200);
  fNTPCIso = 0;
  fTPCIso = new TClonesArray("COMET::IReconTrackerModule::ITPCTrack", 200);
  fNTPCUnused = 0;
  fTPCUnused = new TClonesArray("COMET::IReconTrackerModule::IUnusedHit", 200);
  fNFGDUnused = 0;
  fFGDUnused = new TClonesArray("COMET::IReconTrackerModule::IUnusedHit", 200);
  fNTPCExtrapolation = 0;
  fTPCExtrapolation = new TClonesArray("COMET::IReconTrackerModule::ITPCTrack", 200);
}

//*************************************************************
COMET::IReconTrackerModule::~IReconTrackerModule() {}
//*************************************************************

//*************************************************************
Bool_t COMET::IReconTrackerModule::ProcessFirstEvent(COMET::ICOMETEvent& event) {
//*************************************************************

  fMaxDrift        = COMET::IGeomInfo::Get().TPC().GetMaxDriftDistance(0,0);
  const TGeoNode *cath = COMET::IGeomInfo::Get().TPC().GetCathode(0);
  const TGeoVolume *cathvol =	cath->GetVolume();
  TGeoShape       *cathshape =	cathvol->GetShape();
  TGeoBBox*       box = dynamic_cast<TGeoBBox*>(cathshape);
  fCathodeOffset =  box->GetDX(); 

  return true;
}

//*************************************************************
void COMET::IReconTrackerModule::InitializeModule() {

  fDriftVelocity   = COMET::IOARuntimeParameters::Get().GetParameterD("ReconTPC.TPCgas.DriftSpeed");
 
}
//*************************************************************

//*************************************************************
void COMET::IReconTrackerModule::InitializeBranches() {
  //*************************************************************
  SetSplitLevel(1);
  
  
  // The number of objects that have been saved.
  fOutputTree->Branch("NVertices",
          &fNVertices, "NVertices/I",
          fBufferSize);
  fOutputTree->Branch("NTracks",
          &fNTracks, "NTracks/I",
          fBufferSize);
  fOutputTree->Branch("NFGDOther",
          &fNFGDOther, "NFGDOther/I",
          fBufferSize);
  fOutputTree->Branch("NTPCOther",
          &fNTPCOther, "NTPCOther/I",
          fBufferSize);
  fOutputTree->Branch("NTPCIso",
          &fNTPCIso, "NTPCIso/I",
          fBufferSize);
  fOutputTree->Branch("NFGDUnused",
          &fNFGDUnused, "NFGDUnused/I",
          fBufferSize);
  fOutputTree->Branch("NTPCUnused",
          &fNTPCUnused, "NTPCUnused/I",
          fBufferSize);
  fOutputTree->Branch("NTPCExtrapolation",
          &fNTPCExtrapolation, "NTPCExtrapolation/I",
          fBufferSize);

  // A branch for the clone array of Tracker object information.
  fOutputTree->Branch("Vertices",
          &fVertices,
          fBufferSize, fSplitLevel);
  fOutputTree->Branch("Tracks",
          &fTracks,
          fBufferSize, fSplitLevel);
  fOutputTree->Branch("FGDOther",
          &fFGDOther,
          fBufferSize, fSplitLevel);
  fOutputTree->Branch("TPCOther",
          &fTPCOther,
          fBufferSize, fSplitLevel);
  fOutputTree->Branch("TPCIso",
          &fTPCIso,
          fBufferSize, fSplitLevel);
  fOutputTree->Branch("TPCUnused",
          &fTPCUnused,
          fBufferSize, fSplitLevel);
  fOutputTree->Branch("FGDUnused",
          &fFGDUnused,
          fBufferSize, fSplitLevel);
  fOutputTree->Branch("TPCExtrapolation",
          &fTPCExtrapolation,
          fBufferSize, fSplitLevel);
  return;
}


//*************************************************************
void COMET::IReconTrackerModule::FillVertex(COMET::IHandle< COMET::IReconVertex > in, TClonesArray * out, int idx ){
//*************************************************************

  ITrackerVertex * vertex = new ( (*out)[idx] ) ITrackerVertex;
  
  vertex->AlgorithmName = in->GetAlgorithmName();
  vertex->Status        = in->CheckStatus( in->kSuccess );
  vertex->Quality       = in->GetQuality();
  vertex->NDOF          = in->GetNDOF();
  
  COMET::IHandle< COMET::IVertexState > state = in->GetState();
  if ( state ){
    vertex->Position = state->GetPosition();
    vertex->Variance = state->GetPositionVariance();
  } else {
    vertex->Position = TLorentzVector( 0., 0., 0., 0. );
    vertex->Variance = TLorentzVector( 0., 0., 0., 0. );
  }
  return;
}


//*************************************************************
void COMET::IReconTrackerModule::FillFGD(COMET::IHandle< COMET::IReconTrack > in, TClonesArray * out, int idx ){
//*************************************************************
  
  COMET::IReconBase::Status dets[2] ={COMET::IReconBase::kFGD1, 
           COMET::IReconBase::kFGD2 };

  IFGDTrack * fgd = new ( (*out)[idx] ) IFGDTrack;
  
  for (int i=0; i<2; i++ ){
    if ( in->UsesDetector( dets[i] ) ) fgd->Detector = i+1;
  }
  fgd->UniqueID = in->GetUniqueID();
  fgd->Ndof = in->GetNDOF();
  fgd->Chi2 = in->GetQuality();
  COMET::IHandle< COMET::ITrackState > state = in->GetState();
  if ( state ){
    fgd->Position = state->GetPosition();
    fgd->Direction = state->GetDirection();
    fgd->EDeposit = state->GetEDeposit();
  } else {
    // set bogus values
    fgd->Position = TLorentzVector(0.,0.,0.,0.);
    fgd->Direction = TVector3( 0.,0.,0.);
    fgd->EDeposit = 0.;
  }

  return;
}


//*************************************************************
void COMET::IReconTrackerModule::FillTPC(COMET::IHandle< COMET::IReconPID > in, TClonesArray * out, int idx ){
//*************************************************************
  
  COMET::IReconBase::Status dets[3] ={COMET::IReconBase::kTPC1, 
           COMET::IReconBase::kTPC2,
           COMET::IReconBase::kTPC3 };

  ITPCTrack * tpc = new ( (*out)[idx] ) ITPCTrack;

  for (int i=0; i<3; i++ ){
    if ( in->UsesDetector( dets[i] ) ) tpc->Detector = i+1;
  }
  tpc->Ndof = in->GetNDOF();
  tpc->Chi2 = in->GetQuality();
  tpc->NHits = 0;
  tpc->UniqueID = in->GetUniqueID();
  if( in->GetHits() ) tpc->NHits = GetNumberOfHits( in );

  COMET::IHandle< COMET::IPIDState > state = in->GetState();
  if ( state ){
    tpc->Momentum = state->GetMomentum();
    tpc->MomentumError = sqrt( fabs(state->GetMomentumVariance()) );
    tpc->Charge = state->GetCharge();
    tpc->Position = state->GetPosition();
    tpc->Direction = state->GetDirection();
    tpc->DirectionVariance = state->GetDirectionVariance();
  } else {
    // set bogus values
    tpc->Momentum = 0.;
    tpc->MomentumError = 0.;
    tpc->Charge = 0.;
    tpc->Position = TLorentzVector(0.,0.,0.,0.);
    tpc->Direction = TVector3( 0.,0.,0.);
    tpc->PositionVariance = TLorentzVector(0.,0.,0.,0.);
    tpc->DirectionVariance = TVector3( 0.,0.,0.);
  }

  if( in->Get<COMET::IRealDatum>("tpcPid_NTrun") )
    tpc->NTrun         = in->Get<COMET::IRealDatum>("tpcPid_NTrun")->GetValue();
  else
    tpc->NTrun         = -0xABCDEF;  

  if( in->Get<COMET::IRealDatum>("tpcPid_Ccorr") )
    tpc->Ccorr         = in->Get<COMET::IRealDatum>("tpcPid_Ccorr")->GetValue();
  else
    tpc->Ccorr         = -0xABCDEF;
  
  if ( in->Get<COMET::IRealDatum>("tpcPid_PullEle") )
    tpc->PullEle       = in->Get<COMET::IRealDatum>("tpcPid_PullEle")->GetValue();
  else
    tpc->PullEle       = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_PullMuon") )
    tpc->PullMuon      = in->Get<COMET::IRealDatum>("tpcPid_PullMuon")->GetValue();
  else
    tpc->PullMuon      = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_PullPion") )
    tpc->PullPion      = in->Get<COMET::IRealDatum>("tpcPid_PullPion")->GetValue();
  else
    tpc->PullPion      = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_PullKaon") )
    tpc->PullKaon      = in->Get<COMET::IRealDatum>("tpcPid_PullKaon")->GetValue();
  else
    tpc->PullKaon      = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_PullProton") )
    tpc->PullProton    = in->Get<COMET::IRealDatum>("tpcPid_PullProton")->GetValue();
  else
    tpc->PullProton    = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_dEdxexpEle") )
    tpc->dEdxexpEle    = in->Get<COMET::IRealDatum>("tpcPid_dEdxexpEle")->GetValue();
  else
    tpc->dEdxexpEle    = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_dEdxexpMuon") )
    tpc->dEdxexpMuon   = in->Get<COMET::IRealDatum>("tpcPid_dEdxexpMuon")->GetValue();
  else
    tpc->dEdxexpMuon   = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_dEdxexpPion") )
    tpc->dEdxexpPion   = in->Get<COMET::IRealDatum>("tpcPid_dEdxexpPion")->GetValue();
  else
    tpc->dEdxexpPion   = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_dEdxexpKaon") )
    tpc->dEdxexpKaon   = in->Get<COMET::IRealDatum>("tpcPid_dEdxexpKaon")->GetValue();
  else
    tpc->dEdxexpKaon   = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_dEdxexpProton") )
    tpc->dEdxexpProton = in->Get<COMET::IRealDatum>("tpcPid_dEdxexpProton")->GetValue();
  else
    tpc->dEdxexpProton = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_SigmaEle") )
    tpc->SigmaEle      = in->Get<COMET::IRealDatum>("tpcPid_SigmaEle")->GetValue();
  else
    tpc->SigmaEle      = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_SigmaMuon") )
    tpc->SigmaMuon     = in->Get<COMET::IRealDatum>("tpcPid_SigmaMuon")->GetValue();
  else
    tpc->SigmaMuon     = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_SigmaPion") )
    tpc->SigmaPion     = in->Get<COMET::IRealDatum>("tpcPid_SigmaPion")->GetValue();
  else
    tpc->SigmaPion     = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_SigmaKaon") )
    tpc->SigmaKaon     = in->Get<COMET::IRealDatum>("tpcPid_SigmaKaon")->GetValue();
  else
    tpc->SigmaKaon     = -0xABCDEF;

  if ( in->Get<COMET::IRealDatum>("tpcPid_SigmaProton") )
    tpc->SigmaProton   = in->Get<COMET::IRealDatum>("tpcPid_SigmaProton")->GetValue();
  else
    tpc->SigmaProton   = -0xABCDEF;

  // now look for constituent track that made this PID
  tpc->NConstituents = 0;
  tpc->Sigma0 = 0;
  tpc->Sigma1 = 0;
  tpc->Sigma2 = 0;
  tpc->MeanDrift = 0;
  tpc->TrDirection = TVector3(0,0,0);
  tpc->TrDirectionVar = TVector3(0,0,0);
  tpc->TrCurvature = 0.0;
  tpc->TrCurvatureVar = 0.0;

  // check if this is a composite object... if so fill constituents
  if ( in->GetConstituents() ) {
    COMET::IHandle<COMET::IReconTrack> track = GetConstituentTrack( in, tpc->NConstituents );
    if (track){
      //track->ls(); 
      COMET::IHandle<COMET::ITrackState> astate = track->GetState();
      if (astate){
  tpc->TrDirection    = astate->GetDirection();
  tpc->TrDirectionVar.SetX(  astate->GetCovarianceValue( astate->GetDirectionIndex(),  astate->GetDirectionIndex() ) );
  tpc->TrCurvature    = astate->GetCurvature();
  tpc->TrCurvatureVar = astate->GetCovarianceValue( astate->GetCurvatureIndex(),  astate->GetCurvatureIndex() );
      }
      if ( track->Get<COMET::IRealDatum>("Sigma0") )
  tpc->Sigma0 = track->Get<COMET::IRealDatum>("Sigma0")->GetValue();
      if ( track->Get<COMET::IRealDatum>("Sigma1") )
  tpc->Sigma1 = track->Get<COMET::IRealDatum>("Sigma1")->GetValue();
      if ( track->Get<COMET::IRealDatum>("Sigma2") )
  tpc->Sigma2 = track->Get<COMET::IRealDatum>("Sigma2")->GetValue();
      if ( track->Get<COMET::IRealDatum>("MeanDrift") )
  tpc->MeanDrift = track->Get<COMET::IRealDatum>("MeanDrift")->GetValue();
    }
  }
  tpc->HasExtrapolation = false;
}

// Utility functions from Claudio's code
void getTangents(COMET::IHandle<COMET::ITrackState> state,double &tx,double &etx,double &ty,double &ety){

  tx = state->GetDirection()[0]/state->GetDirection()[2]; 
  ty = state->GetDirection()[1]/state->GetDirection()[2]; 

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

  for(int i = 0; i < 3; i++ ) 
    for(int j = 0; j < 3; j++ ) {
      Covariance[i][j] = state->GetCovarianceValue(state->GetDirectionIndex()+i,state->GetDirectionIndex()+j);
    }

  Jacobian[0][0] = 1./state->GetDirection()[2]; 
  Jacobian[0][1] = 0.;
  Jacobian[0][2] = -tx/state->GetDirection()[2];

  Jacobian[1][0] = 0.;
  Jacobian[1][1] = 1./state->GetDirection()[2]; 
  Jacobian[1][2] = -ty/state->GetDirection()[2];
  
  Cov = Jacobian*Covariance*Jacobian.T();

  etx = TMath::Sqrt(Cov[0][0]);
  ety = TMath::Sqrt(Cov[1][1]);
}

//=========================================================//
bool insideFGD1(double zx, double zy, double xx, double yy)
{
  return ((zx>140)&&(zy>140)&&(zx<500)&&(zy<500)&&(fabs(yy)<800)&&(fabs(xx)<800));
}
//=========================================================//
bool insideFGD2(double zx, double zy, double xx, double yy)
{
  return ((zx>1530)&&(zy>1530)&&(zx<1900)&&(zy<1900)&&(fabs(yy)<800)&&(fabs(xx)<800));
}
//=========================================================//
bool borderFGD1(TVector3 pos)
{
  bool ret=true;
  if((pos.Y() < 850) && (pos.Y() > -800) && (pos.X()<850))
    ret = false;
  if(pos.Z()>500)
    ret = false;
  if(fabs(pos.Z()-120.85)<0.0001){
    ret=true;
  }
  return ret;
}
//=========================================================//
bool borderFGD2(TVector3 pos){
  bool ret=true;
  if((pos.Y()<850)&&(pos.Y()>-800)&&(fabs(pos.X())<850))
    {ret=false;}
  if(pos.Z()<500)
    ret=false;
  if(fabs(pos.Z()-1478.85)<0.0001)
    ret=true;
  return ret;
}

// This is the distance computation function used in the extrapolation analysis
double computeDistFGD(TVector3 f, TVector3 in, double tx, double ty, double curv, bool isX)
{
  // reject hits far away, tracks with zero curvature, and hits downstream from the track
  if ((fabs(curv)<1.e-20) || (fabs(in.Z()-f.Z())>600) || (f.Z()>in.Z()))
    return -9999.0;
  
  double r1 = fabs(1./curv);
  double phi1 = atan(ty);
  double charge1 = 1.;
  if (curv>0) charge1=-1.;
  double zc1 = in.Z()-charge1*r1*sin(phi1);
  double yc1 = in.Y()+charge1*r1*cos(phi1);
  double delta =  pow(r1,2) - pow((zc1-f.Z()),2) ;
  if (delta<0) 
    return -9999.0; 
  double ySol1 = yc1 + sqrt(delta);
  double ySol2 = yc1 - sqrt(delta);
  double ySol=ySol1;
  if (fabs(ySol1-in.Y()) > fabs(ySol2-in.Y()) ) ySol=ySol2;
  // check with linear extrapo

  // yInts currentl not being used, so commented out to remove compiler warnings
  //double yIntS = in.Y()+ty*(f.Z()- in.Z());

  if (!isX) {
    return f.Y()-ySol;
  } else {
    return f.X() - in.X() - tx*(f.Z()-in.Z());
  }
}

//*************************************************************
bool COMET::IReconTrackerModule::ExtrapolateTPC(COMET::IHandle< COMET::IReconPID > in, TClonesArray * out, int idx, COMET::IHandle<COMET::IHitSelection> fgd )
{
  COMET::IHandle<COMET::IReconTrack> track = in->GetConstituents()->front();
  COMET::IReconNodeContainer &nodes = track->GetNodes();
      
  // Only do the extrapolation matching if there are enough TPC points (36 for now)
  if (nodes.size() < 36)
    return false;
  // ignore TPC1 tracks
  if (in->UsesDetector(COMET::IReconBase::kTPC1))
    return false;

  COMET::IHandle<COMET::ITrackState> frontstate = nodes.front()->GetState();
  if (!frontstate)
    return false;
  if (!fgd)
    return false;

  double txin,tyin,etxin,etyin;
  getTangents(frontstate,txin,etxin,tyin,etyin);
  TVector3 pos = frontstate->GetPosition().Vect();
  double curvature = track->GetCurvature();
  
  // start with bogus numbers
  double ExtrapolatedVertexXX = 9999;
  double ExtrapolatedVertexZY = 9999;
  double ExtrapolatedVertexZX = 9999;
  double ExtrapolatedVertexYY = 9999;
  // std::vector<TVector3*> ExtrapolatedHits;

  // calculate the vertex and find the set of associated FGD hits
  // the distance calculation function will reject hits from the wrong fgd
  // by discarding hits downstream of the track origin and more than 60 cm away in Z
  for (COMET::IHitSelection::iterator hit = fgd->begin();
       hit != fgd->end(); hit++) {
    TVector3 fgdPos = (*hit)->GetPosition();
  
    double dist = computeDistFGD(fgdPos, pos, txin, tyin, curvature,(*hit)->IsXHit());
    if (fabs(dist) >= 30) // we reject this hit (30 mm is the threshold)
      continue;

    // ExtrapolatedHits.push_back(new TVector3(fgdPos));

    if ((*hit)->IsYHit() && fgdPos.Z() < ExtrapolatedVertexZY) {  
      ExtrapolatedVertexZY = fgdPos.Z();
      ExtrapolatedVertexYY = fgdPos.Y();
    }//Y coord
  
    if ((*hit)->IsXHit() && fgdPos.Z() < ExtrapolatedVertexZX){
      ExtrapolatedVertexZX = fgdPos.Z();
      ExtrapolatedVertexXX = fgdPos.X();
    }// X coord
  }// loop on FGD hits
  
  if (ExtrapolatedVertexZX == 9999 || ExtrapolatedVertexZY == 9999)
    return false;

  bool EnteringFGD = false;
  // Now check whether this qualifies as an 'entering' track
  if (insideFGD1(ExtrapolatedVertexZX,ExtrapolatedVertexZY,ExtrapolatedVertexXX,ExtrapolatedVertexYY)) {
    for (COMET::IHitSelection::iterator hit = fgd->begin();
   hit != fgd->end(); hit++) {
      TVector3 fgdPos = (*hit)->GetPosition();
  
      if (!borderFGD1(fgdPos))
  continue;

      double dist = computeDistFGD(fgdPos,pos,txin,tyin,curvature,(*hit)->IsXHit());
      if (fabs(dist) < 150) { // cutoff here is 150 mm
  EnteringFGD = true;
  break;
      }
    }
  }
    
  if(insideFGD2(ExtrapolatedVertexZX,ExtrapolatedVertexZY,ExtrapolatedVertexXX,ExtrapolatedVertexYY)){
    for (COMET::IHitSelection::iterator hit = fgd->begin();
   hit != fgd->end(); hit++) {
      TVector3 fgdPos = (*hit)->GetPosition();

      if (!borderFGD2(fgdPos))
  continue;

      double dist = computeDistFGD(fgdPos,pos,txin,tyin,curvature,(*hit)->IsXHit());
      if (fabs(dist) < 150) {
  EnteringFGD = true;
  break;
      }
    }
  }

  // now we fill the rest of the info for the object
  FillTPC(in, out, idx);
  ITPCTrack* tpc = dynamic_cast<ITPCTrack*>((*out)[idx]);
  if (!tpc)
    return false;

  tpc->EnteringFGD = EnteringFGD;
  tpc->ExtrapolatedVertexYY = ExtrapolatedVertexYY;
  tpc->ExtrapolatedVertexXX = ExtrapolatedVertexXX;
  tpc->ExtrapolatedVertexZY = ExtrapolatedVertexZY;
  tpc->ExtrapolatedVertexZX = ExtrapolatedVertexZX;
  // std::copy(ExtrapolatedHits.begin(), ExtrapolatedHits.end(), std::back_inserter(tpc->ExtrapolatedHits));
  tpc->HasExtrapolation = true;
  return true;


}
//*************************************************************

/// called recursively to look for bottom level track that built a PID.  In case
/// the PID is made of several PIDs, the one with the track with smallest relative curvature
/// error is returned.
COMET::IHandle< COMET::IReconTrack> COMET::IReconTrackerModule::GetConstituentTrack( COMET::IHandle< COMET::IReconPID > in, int &nconstit ){
  
  COMET::IHandle<COMET::IReconTrack> retval = in;

  if ( in->GetConstituents() ){
    //std::cout<<"In GetConstituentTrack, nconstit="<<nconstit<<std::endl;
    COMET::IReconObjectContainer::const_iterator it;
    for (it=in->GetConstituents()->begin(); it != in->GetConstituents()->end(); it++){
      nconstit++;
      COMET::IHandle<COMET::IReconTrack> track = (*it);
      COMET::IHandle<COMET::IReconPID>   pid = (*it);
      if ( pid ){
  //std::cout<<"  Found pid, dig deeper"<<std::endl;
  COMET::IHandle<COMET::IReconTrack> atrack = GetConstituentTrack( pid, nconstit );
  //if( atrack){
    //std::cout<<"  deeper, found track"<<std::endl;
    //atrack->ls();
  //}
  if ( retval && atrack){
    // compare atrack to retval, save best to retval
    COMET::IHandle<COMET::ITrackState> astate = atrack->GetState();
    COMET::IHandle<COMET::ITrackState> rstate = retval->GetState();
    if ( astate && !rstate ) retval = atrack;
    if ( rstate && astate ){
      double amom, amomerr;
      double rmom, rmomerr;
      GetMomentum( astate, amom, amomerr );
      GetMomentum( rstate, rmom, rmomerr );
      if ( amom > 0.0 && rmom <= 0.0 ) retval = atrack;
      if ( amom > 0.0 && rmom > 0.0 && amomerr/amom < rmomerr/rmom ) retval = atrack;
    }
  } else if ( atrack ){
    retval = atrack;
  }
      } else if ( track ){
  //std::cout<<"  Found track"<<std::endl;
  //track->ls();
  if (retval){
    // compare track to retval
    // compare atrack to retval, save best to retval
    COMET::IHandle<COMET::ITrackState> astate = track->GetState();
    COMET::IHandle<COMET::ITrackState> rstate = retval->GetState();
    if ( astate && !rstate ) retval = track;
    if ( rstate && astate ){
      double amom, amomerr;
      double rmom, rmomerr;
      GetMomentum( astate, amom, amomerr );
      GetMomentum( rstate, rmom, rmomerr );
      if ( amom > 0.0 && rmom <= 0.0 ) retval = track;
      if ( amom > 0.0 && rmom > 0.0 && amomerr/amom < rmomerr/rmom ) retval = track;
    }
  } else {
    retval = track;
  }
      }
    }
  }
  
  return retval;
}



//**************************************************************
COMET::IHandle<COMET::IHitSelection> COMET::IReconTrackerModule::GetHitsFromComboHits(COMET::IHandle<COMET::IHitSelection> tpc) {
  //**************************************************************

  if (!tpc) return COMET::IHandle<COMET::IHitSelection>();

  COMET::IHandle<COMET::IHitSelection> hits(new COMET::IHitSelection());

  for(COMET::IHitSelection::const_iterator hsel = tpc->begin(); hsel != tpc->end(); hsel++){
    COMET::IHandle<COMET::IComboHit> combo(*hsel);
    if (combo){
      for (COMET::IHitSelection::const_iterator cHit = combo->GetHits().begin();
     cHit != combo->GetHits().end();cHit++){
  if (*cHit)
    hits->push_back(*cHit);
      }
    }
  }

  return hits;
  
}



//*************************************************************
void COMET::IReconTrackerModule::FillOther(COMET::IHandle< COMET::IReconTrack > in, TClonesArray * out, int idx ){
//*************************************************************
  
  COMET::IReconBase::Status dets[5] ={COMET::IReconBase::kTPC1, 
           COMET::IReconBase::kTPC2,
           COMET::IReconBase::kTPC3,
           COMET::IReconBase::kFGD1,
           COMET::IReconBase::kFGD2};

  ITrackOther* other = new ( (*out)[idx] ) ITrackOther;

  other->AlgorithmName = in->GetAlgorithmName();
  other->Detector = -1;
  for (int i=0; i<5; i++ ){
    if ( in->UsesDetector( dets[i] ) ) other->Detector = i+1;
  }
  COMET::IHandle< COMET::IHitSelection > hits = GetHitsFromComboHits( in->GetHits() );
  // get track t0:
  double t0 = in->GetPosition().T();

  other->EDeposit = 0.0;

  //Initialise mint and maxt to remove compiler warnings

  TVector3 minpos(0,0,999999);  double mint = 0;
  TVector3 maxpos(0,0,-999999);  double maxt = 0;
  for ( COMET::IHitSelection::const_iterator hsel = hits->begin(); hsel != hits->end(); hsel++){
    IUnusedHit * hit = new ( (*(other->Hits))[other->NHits++] ) IUnusedHit;
    //    hit->TotalCharge = (*hsel)->GetCharge();
    hit->Position    = (*hsel)->GetPosition();
    hit->Variance    = (*hsel)->GetUncertainty();
    //hit->Time        = (*hsel)->GetTime();
    hit->TimeUnc     = (*hsel)->GetTimeUncertainty();

    double peakTime=0;
    double maxCharge = 0; 
    
    // Find the pick charge and Time of the Wave form
    COMET::IHandle<COMET::IMultiHit> mh = *hsel;
    std::vector< COMET::IHandle< COMET::ISingleHit > >::const_iterator mhit;
    for( mhit = mh->begin(); mhit !=  mh->end(); mhit++ ) {
      if( maxCharge < (*mhit)->GetCharge() ) {
  maxCharge = (*mhit)->GetCharge(); 
  peakTime = (*mhit)->GetTime();
      }
    }
    
    hit->TotalCharge = maxCharge;

    COMET::IGeometryId geomId = (*hsel)->GetGeomId();
    double t = peakTime - t0;
    double x = (fMaxDrift + fCathodeOffset - t*fDriftVelocity)*COMET::IGeomInfo::Get().TPC().GetDriftSense(geomId);

    hit->Position[0] = x;
    hit->Time = t;

    other->EDeposit += hit->TotalCharge;
    if ( minpos[2] > hit->Position[2] ){
      minpos = hit->Position; 
      mint   = hit->Time;
    }
    if ( maxpos[2] < hit->Position[2] ){
      maxpos = hit->Position; 
      maxt   = hit->Time;
    }
  }

  other->FrontPosition = TLorentzVector( minpos[0], minpos[1], minpos[2], mint );
  other->BackPosition = TLorentzVector( maxpos[0], maxpos[1], maxpos[2], maxt );

  return;
}


    
//*************************************************************
// recursively called to fill all constituents, constituents of constituents, and so on
// if new constituent is created return true, otherwise return false
bool COMET::IReconTrackerModule::FillConstituentPid( ITrackerResult * result, COMET::IHandle< COMET::IReconPID > in, int &idx ){
//*************************************************************

  // First check if this is single detector pid (again only check TPC, since FGD is track)
  std::string fDetName = std::string( in->ConvertDetector() );
  if ( fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)" ){
    FillTPC( in , result->TPC, result->NTPCs );
    result->NTPCs++;  
    return false;
  }

  // otherwise it is indeed a composite object so create new constituent object
  ITrackerConstituent * constituent = new ( (*(result->Constituents))[result->NTotalConstituents++] ) ITrackerConstituent;
  idx = result->NTotalConstituents - 1 ;
  constituent->AlgorithmName = in->GetAlgorithmName();
  constituent->Detectors     = GetDetectorNumber( in );
  constituent->Status        = in->CheckStatus( in->kSuccess );
  constituent->NDOF          = in->GetNDOF();
  constituent->Chi2          = in->GetQuality();
  if (constituent->NDOF>0 ) constituent->Quality = TMath::Prob( constituent->Chi2, constituent->NDOF);
  else result->Quality = -1.0;
    
  // fill constituent state information
  COMET::IHandle<COMET::IPIDState> state = in->GetState();
  if( state ){
    constituent->Position      = state->GetPosition();
    constituent->Variance      = state->GetPositionVariance();
    constituent->Direction     = state->GetDirection();
    constituent->isForward     = ( constituent->Direction.Z() > 0 );
    constituent->Momentum      = state->GetMomentum();
    constituent->MomentumError = sqrt( fabs( state->GetMomentumVariance() ) );
    constituent->Charge        = state->GetCharge();
  } else {
    constituent->Position      = TLorentzVector( 0., 0., 0., 0. );	
    constituent->Variance      = TLorentzVector( 0., 0., 0., 0. );	
    constituent->Direction     = TVector3( 0., 0., 0. );
    constituent->isForward     = 1;
    constituent->Momentum      = 0.0;
    constituent->MomentumError = 0.0;
    constituent->Charge        = 0.0;
  }
  constituent->EDeposit      = 0.0;  //No edeposit variable for pid recon objects


  // Get state at front and back of track
  constituent->NNodes = in->GetNodes().size();
  if (constituent->NNodes>0){
    COMET::IHandle<COMET::IPIDState> firstState = in->GetNodes()[0]->GetState();
    COMET::IHandle<COMET::IPIDState> lastState = in->GetNodes()[in->GetNodes().size()-1]->GetState();
    
    if ( firstState ){
      constituent->FrontPosition  = firstState->GetPosition();
      constituent->FrontVariance  = firstState->GetPositionVariance();
      constituent->FrontDirection = firstState->GetDirection();
      constituent->FrontMomentum  = firstState->GetMomentum();
    } else {
      constituent->FrontPosition  = TLorentzVector( 0., 0., 0., 0. );	
      constituent->FrontVariance  = TLorentzVector( 0., 0., 0., 0. );
      constituent->FrontDirection = TVector3( 0., 0., 0. );
      constituent->FrontMomentum  = 0.;
    }
    
    if ( lastState ){
      constituent->BackPosition  = lastState->GetPosition();
      constituent->BackVariance  = lastState->GetPositionVariance();
      constituent->BackDirection = lastState->GetDirection();
      constituent->BackMomentum  = lastState->GetMomentum();
    } else {
      constituent->BackPosition  = TLorentzVector( 0., 0., 0., 0. );	
      constituent->BackVariance  = TLorentzVector( 0., 0., 0., 0. );
      constituent->BackDirection = TVector3( 0., 0., 0. );
      constituent->BackMomentum  = 0.;
    }
  }
  
  // get number of hits if available
  constituent->NHits = 0;
  if ( in->GetHits() ) constituent->NHits = GetNumberOfHits( in );
  
  //Check if this is a single detector object?  If it is, then just save it, and don't worry about 
  // any constituents inside of single detector
  constituent->NConstituents = 0;
  // check if this is a composite object... if so fill constituents
  if ( in->GetConstituents() ){
    COMET::IReconObjectContainer::const_iterator it;
    for (it=in->GetConstituents()->begin(); it != in->GetConstituents()->end(); it++){
      COMET::IHandle<COMET::IReconTrack> track = (*it);
      COMET::IHandle<COMET::IReconPID>   pid = (*it);
      if ( pid ){
  if ( constituent->NConstituents >= 2 ) {
    //std::cout<<"IReconTrackerModule::FillConstituentPID>pid nconstituents already 2, not adding more"<<std::endl;
  } else {
    bool addedconstituent = FillConstituentPid( result, pid, constituent->ConstitIdx[ constituent->NConstituents ] );
    if (addedconstituent) constituent->NConstituents++;
  }
      } else if ( track ){
  if ( constituent->NConstituents >= 2 ) {
    // std::cout<<"IReconTrackerModule::FillConstituentPID>track nconstituents already 2, not adding more"<<std::endl;
  } else {
    bool addedconstituent = FillConstituentTrack( result, track, constituent->ConstitIdx[ constituent->NConstituents ] );
    if (addedconstituent) constituent->NConstituents++;
  }
      }
    }
  }
  
  return true;
}


    
//*************************************************************
// recursively called to fill all constituents, constituents of constituents, and so on
// if new constituent is created return true, otherwise return false
bool COMET::IReconTrackerModule::FillConstituentTrack( ITrackerResult * result, COMET::IHandle< COMET::IReconTrack > in, int& idx ){
//*************************************************************

  // First check if this is single detector
  std::string fDetName = std::string( in->ConvertDetector() );
  if ( fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)" ){
    FillTPC( in , result->TPC, result->NTPCs );
    result->NTPCs++;  
    return false;
  }
  if ( fDetName == "(FGD1)" || fDetName == "(FGD2)" ){
    FillFGD( in , result->FGD, result->NFGDs );
    result->NFGDs++;  
    return false;
  }

  // otherwise it is indeed a composite object so create new constituent object
  ITrackerConstituent * constituent = new ( (*(result->Constituents))[result->NTotalConstituents++] ) ITrackerConstituent;
  idx = result->NTotalConstituents - 1 ;
  constituent->AlgorithmName = in->GetAlgorithmName();
  constituent->Detectors     = GetDetectorNumber( in );
  constituent->Status        = in->CheckStatus( in->kSuccess );
  constituent->NDOF          = in->GetNDOF();
  constituent->Chi2          = in->GetQuality();
  if (constituent->NDOF>0 ) constituent->Quality = TMath::Prob( constituent->Chi2, constituent->NDOF);
  else result->Quality = -1.0;
    
  // fill constituent state information
  COMET::IHandle<COMET::ITrackState> state = in->GetState();
  if( state ){
    constituent->Position      = state->GetPosition();
    constituent->Variance      = state->GetPositionVariance();
    constituent->Direction     = state->GetDirection();
    constituent->DirectionVariance = state->GetDirectionVariance();
    constituent->isForward     = ( constituent->Direction.Z() > 0 );
    GetMomentum( state, constituent->Momentum, constituent->MomentumError );
    constituent->Charge        = GetCharge( state );
    constituent->EDeposit      = state->GetEDeposit(); 
  } else {
    constituent->Position      = TLorentzVector( 0., 0., 0., 0. );	
    constituent->Variance      = TLorentzVector( 0., 0., 0., 0. );	
    constituent->Direction     = TVector3( 0., 0., 0. );
    constituent->DirectionVariance = TVector3( 0., 0., 0. );
    constituent->isForward     = 1;
    constituent->Momentum      = 0.0;
    constituent->MomentumError = 0.0;
    constituent->Charge        = 0.0;
  }

  // Get state at front and back of track
  constituent->NNodes = in->GetNodes().size();
  if (constituent->NNodes>0){
    COMET::IHandle<COMET::ITrackState> firstState = in->GetNodes()[0]->GetState();
    COMET::IHandle<COMET::ITrackState> lastState = in->GetNodes()[in->GetNodes().size()-1]->GetState();
    
    if ( firstState ){
      constituent->FrontPosition  = firstState->GetPosition();
      constituent->FrontVariance  = firstState->GetPositionVariance();
      constituent->FrontDirection = firstState->GetDirection();
      double tmperr;
      GetMomentum( firstState, constituent->FrontMomentum, tmperr );
    } else {
      constituent->FrontPosition  = TLorentzVector( 0., 0., 0., 0. );	
      constituent->FrontVariance  = TLorentzVector( 0., 0., 0., 0. );
      constituent->FrontDirection = TVector3( 0., 0., 0. );
      constituent->FrontMomentum  = 0.;
    }
    
    if ( lastState ){
      constituent->BackPosition  = lastState->GetPosition();
      constituent->BackVariance  = lastState->GetPositionVariance();
      constituent->BackDirection = lastState->GetDirection();
      double tmperr;
      GetMomentum( lastState, constituent->BackMomentum, tmperr );
    } else {
      constituent->BackPosition  = TLorentzVector( 0., 0., 0., 0. );	
      constituent->BackVariance  = TLorentzVector( 0., 0., 0., 0. );
      constituent->BackDirection = TVector3( 0., 0., 0. );
      constituent->BackMomentum  = 0.;
    }
  }
  
  // get number of hits if available
  constituent->NHits = 0;
  if ( in->GetHits() ) constituent->NHits = GetNumberOfHits( in );
  
  //Check if this is a single detector object?  If it is, then just save it, and don't worry about 
  // any constituents inside of single detector
  constituent->NConstituents = 0;
  // check if this is a composite object... if so fill constituents
  if ( in->GetConstituents() ){
    COMET::IReconObjectContainer::const_iterator it;
    for (it=in->GetConstituents()->begin(); it != in->GetConstituents()->end(); it++){
      COMET::IHandle<COMET::IReconTrack> track = (*it);
      COMET::IHandle<COMET::IReconPID>   pid = (*it);
      if ( pid ){
  if ( constituent->NConstituents >= 2 ) {
    //std::cout<<"IReconTrackerModule::FillConstituentTrack>pid nconstituents already 2, not adding more"<<std::endl;
  } else {
    bool addedconstituent = FillConstituentPid( result, pid, constituent->ConstitIdx[ constituent->NConstituents ] );
    if (addedconstituent) constituent->NConstituents++;
  }
      } else if ( track ){
  if ( constituent->NConstituents >= 2 ) {
    //std::cout<<"IReconTrackerModule::FillConstituentPID>track nconstituents already 2, not adding more"<<std::endl;
  } else {	  
    bool addedconstituent = FillConstituentTrack( result, track, constituent->ConstitIdx[ constituent->NConstituents ] );
    if (addedconstituent) constituent->NConstituents++;
  }
      }
    }
  }

  return true;
}


    
//*************************************************************
void COMET::IReconTrackerModule::FillFinalPid( ITrackerResult * result, COMET::IHandle< COMET::IReconPID > reco ){
//*************************************************************

  // Get state at vertex
  COMET::IHandle<COMET::IPIDState> state = reco->GetState();
  if (state){
    result->Position      = state->GetPosition();
    result->Variance      = state->GetPositionVariance();
    result->Direction     = state->GetDirection();
    result->DirectionVariance= state->GetDirectionVariance();
    result->isForward     = ( result->Direction.Z() > 0 );
    result->Momentum      = state->GetMomentum();
    result->MomentumError = sqrt( fabs( state->GetMomentumVariance() ) );
    result->Charge        = state->GetCharge();
  } else {
    result->Position      = TLorentzVector( 0., 0., 0., 0. );	
    result->Variance      = TLorentzVector( 0., 0., 0., 0. );	
    result->Direction     = TVector3( 0., 0., 0. );
    result->DirectionVariance = TVector3( 0., 0., 0. );
    result->isForward     = 1;
    result->Momentum      = 0.;
    result->MomentumError = 0.;
    result->Charge        = 0.;
  }
  result->EDeposit      = 0.;  //No edeposit variable for pid recon objects

  // Get state at front and back of track
  result->NNodes = reco->GetNodes().size();
  if (result->NNodes>0){
    // Loop over nodes
    for (int inode=0; inode<result->NNodes; inode++){
      COMET::IHandle<COMET::IPIDState> iState = reco->GetNodes()[ inode ]->GetState();
      ITrackerNode * node = new ( (*(result->Nodes))[inode] ) ITrackerNode;
      node->EDeposit=0.0;
      if ( iState ){
  node->Charge    = iState->GetCharge();
  node->Position  = iState->GetPosition();
  node->Variance  = iState->GetPositionVariance();
  node->Direction = iState->GetDirection();
  node->DirectionVariance = iState->GetDirectionVariance();
  node->Momentum  = iState->GetMomentum();
  node->MomentumError = sqrt( fabs( iState->GetMomentumVariance() ) );
      } else {
  node->Charge    = 0.0;
  node->Position  = TLorentzVector( 0., 0., 0., 0. );	
  node->Variance  = TLorentzVector( 0., 0., 0., 0. );
  node->Direction = TVector3( 0., 0., 0. );
  node->DirectionVariance = TVector3( 0., 0., 0. );
  node->Momentum  = 0.;
  node->MomentumError = 0.;
      }      
    }
  }
  
  // use recpack to find length of track
  Trajectory traj;
  COMET::IHandle<COMET::IReconBase> object = reco;
  bool ok = COMET::converter().TReconBase_to_Trajectory(object,traj);
  result->Length = 0;
  if (ok){
    double length = 0;
    ok = COMET::rpman().matching_svc().compute_length(traj, length);
    if (ok) result->Length = length;
  }
  
  // get number of hits if available
  result->NHits = 0;
  if ( reco->GetHits() ) result->NHits = GetNumberOfHits( reco );
  
  // get Pid information
  result->Likelihoods.push_back( reco->GetPIDWeight() );
  result->Pids.push_back( reco->GetParticleId() );
  for (COMET::IReconObjectContainer::const_iterator it=reco->GetAlternates().begin(); it!=reco->GetAlternates().end(); it++){
      COMET::IHandle<COMET::IReconPID> alter = *it;
      result->Likelihoods.push_back(alter->GetPIDWeight());
      result->Pids.push_back(alter->GetParticleId());
  }

  // Check if this is a single detector object?  If it is, then just save it, and don't worry about 
  // any constituents inside of single detector
  result->NConstituents = 0;
  std::string fDetName = std::string( reco->ConvertDetector() );
  if ( fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)" ){
    FillTPC( reco , result->TPC, result->NTPCs );
    result->NTPCs++;
  } else {
    // Note that FGD only occurs as IReconTrack, so no need to check that here.
    // Instead this must be a composite object... so fill constituents
    if ( reco->GetConstituents() ){
      COMET::IReconObjectContainer::const_iterator it;
      int iconstit=0;
      for (it=reco->GetConstituents()->begin(); it != reco->GetConstituents()->end(); it++){
  COMET::IHandle<COMET::IReconTrack> track = (*it);
  COMET::IHandle<COMET::IReconPID>   pid = (*it);
  if ( pid ){
    if ( result->NConstituents >= 2 ) {
      //std::cout<<"IReconTrackerModule::FillFinalPid>pid nconstituents already 2, not adding more"<<std::endl;
    } else {
      bool addedconstituent = FillConstituentPid( result, pid, result->ConstitIdx[ result->NConstituents ] );
      if (addedconstituent) result->NConstituents++;
    }
  } else if ( track ){
    if ( result->NConstituents >= 2 ) {
      //std::cout<<"IReconTrackerModule::FillFinalPid>track nconstituents already 2, not adding more"<<std::endl;
    } else {
      bool addedconstituent = FillConstituentTrack( result, track, result->ConstitIdx[ result->NConstituents] );
      if (addedconstituent) result->NConstituents++;
    }
  }
      }
      iconstit++;
    }
  }

  return;
}


//*************************************************************
void COMET::IReconTrackerModule::GetMomentum( COMET::IHandle<COMET::ITrackState> state,double &p,double &ep){
//*************************************************************

  // need to read in magnet field...
  double B =  COMET::IFieldManager::GetFieldValue( state->GetPosition() ).Mag()/unit::tesla;

  // should get speed of light constant ~0.3 from somewhere
  p = 0.3*B/(state->GetCurvature()*TMath::Sqrt(1.-state->GetDirection()[0]*state->GetDirection()[0]));

  CLHEP::HepMatrix  Jacobian(2,2,0);
  CLHEP::HepMatrix  Covariance(2,2,0);
  CLHEP::HepMatrix  Cov(2,2,0);
  
  Covariance[0][0] = state->GetCovarianceValue(state->GetDirectionIndex(),state->GetDirectionIndex());
  Covariance[0][1] = state->GetCovarianceValue(state->GetDirectionIndex(),state->GetCurvatureIndex());

  Covariance[1][0] = state->GetCovarianceValue(state->GetCurvatureIndex(),state->GetDirectionIndex());
  Covariance[1][1] = state->GetCovarianceValue(state->GetCurvatureIndex(),state->GetCurvatureIndex());
  
  Jacobian[0][0] = 1.;
  Jacobian[0][1] = 0.;

  Jacobian[1][0] = p/(1.-state->GetDirection()[0]*state->GetDirection()[0])*TMath::Abs(state->GetDirection()[0]);
  Jacobian[1][1] = -p/state->GetCurvature(); 

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

  ep = TMath::Sqrt(Cov[1][1]);

  p = TMath::Abs(p);
  
  return ; 
}



//*************************************************************
double COMET::IReconTrackerModule::GetCharge( COMET::IHandle<COMET::ITrackState> state){
//*************************************************************
  double q = ( state->GetCurvature()>0.0? -1.0:1.0 );
  if ( state->GetDirection()[2] < 0.0 ) q*=-1.0;
  return q;
}


//*************************************************************
void COMET::IReconTrackerModule::FillFinalTrack( ITrackerResult * result, COMET::IHandle< COMET::IReconTrack > reco ){
//*************************************************************
  
  // Get state at vertex
  COMET::IHandle<COMET::ITrackState> state = reco->GetState();
  if (state){
    result->Position      = state->GetPosition();
    result->Variance      = state->GetPositionVariance();
    result->Direction     = state->GetDirection();
    result->DirectionVariance = state->GetDirectionVariance();
    result->isForward     = ( result->Direction.Z() > 0 );
    result->EDeposit      = state->GetEDeposit();
    // need conversion from curvature to momentum
    GetMomentum( state, result->Momentum, result->MomentumError );
    // need conversion from curvature, and direction to charge
    result->Charge        = GetCharge( state );
  } else {
    result->Position      = TLorentzVector( 0., 0., 0., 0. );	
    result->Variance      = TLorentzVector( 0., 0., 0., 0. );	
    result->Direction     = TVector3( 0., 0., 0. );
    result->DirectionVariance = TVector3( 0., 0., 0. );
    result->isForward     = 1;
    result->Momentum      = 0.;
    result->MomentumError = 0.;
    result->Charge        = 0.;
  }

  // Get state at front and back of track
  result->NNodes = reco->GetNodes().size();
  if (result->NNodes>0){
    // Loop over nodes
    for (int inode=0; inode<result->NNodes; inode++){
      COMET::IHandle<COMET::ITrackState> iState = reco->GetNodes()[ inode ]->GetState();
      ITrackerNode * node = new ( (*(result->Nodes))[inode] ) ITrackerNode;
      if ( iState ){
  node->EDeposit  = iState->GetEDeposit();
  node->Position  = iState->GetPosition();
  node->Variance  = iState->GetPositionVariance();
  node->Direction = iState->GetDirection();
  node->DirectionVariance = iState->GetDirectionVariance();
  // need conversion from curvature, and direction to charge
  GetMomentum( iState, node->Momentum, node->MomentumError );
  node->Charge     = GetCharge( iState );
      } else {
  node->EDeposit  = 0.0;
  node->Charge    = 0.0;
  node->Position  = TLorentzVector( 0., 0., 0., 0. );	
  node->Variance  = TLorentzVector( 0., 0., 0., 0. );
  node->Direction = TVector3( 0., 0., 0. );
  node->DirectionVariance = TVector3( 0., 0., 0. );
  node->Momentum  = 0.;
  node->MomentumError = 0.;
      }      
    }
  }
  

  // use recpack to find length of track
  Trajectory traj;
  COMET::IHandle<COMET::IReconBase> object = reco;
  bool ok = COMET::converter().TReconBase_to_Trajectory(object,traj);
  result->Length = 0;
  if (ok){
    double length = 0;
    ok = COMET::rpman().matching_svc().compute_length(traj, length);
    if (ok) result->Length = length;
  }

  // get number of hits if available
  result->NHits = 0;
  if ( reco->GetHits() ) result->NHits = GetNumberOfHits( reco );
  
  //Check if this is a single detector object?  If it is, then just save it, and don't worry about 
  // any constituents inside of single detector
  result->NConstituents = 0;
  std::string fDetName = std::string( reco->ConvertDetector() );
  if ( fDetName == "(TPC1)" || fDetName == "(TPC2)" || fDetName == "(TPC3)" ){
    // I don't think I should ever get here, since I should have found IReconPID for the TPC instead
    FillTPC( reco , result->TPC, result->NTPCs );
    result->NTPCs++;
  } else if ( fDetName == "(FGD1)" || fDetName == "(FGD2)" ){
    FillFGD( reco, result->FGD, result->NFGDs );
    result->NFGDs++;
  } else {
     // Instead this must be a composite object... so fill constituents
    if ( reco->GetConstituents() ){
      COMET::IReconObjectContainer::const_iterator it;
      for (it=reco->GetConstituents()->begin(); it != reco->GetConstituents()->end(); it++){
  COMET::IHandle<COMET::IReconTrack> track = (*it);
  COMET::IHandle<COMET::IReconPID>   pid = (*it);
  if ( pid ){
    if ( result->NConstituents >= 2 ) {
      //std::cout<<"IReconTrackerModule::FillFinalTrack>pid nconstituents already 2, not adding more"<<std::endl;
    } else {
      bool addedconstituent = FillConstituentPid( result, pid, result->ConstitIdx[ result->NConstituents ] );
      if (addedconstituent) result->NConstituents++;
    }
  } else if ( track ){
    if ( result->NConstituents >= 2 ) {
      //std::cout<<"IReconTrackerModule::FillFinalPid>track nconstituents already 2, not adding more"<<std::endl;
    } else {
      bool addedconstituent = FillConstituentTrack( result, track, result->ConstitIdx[ result->NConstituents ] );
      if (addedconstituent) result->NConstituents++;
    }
  }
      }
    }
  }

  return;
}


//*****************************************************
int COMET::IReconTrackerModule::GetDetectorNumber(COMET::IHandle<COMET::IReconBase> object){
//*****************************************************

  int det=0;

  if (object->UsesDetector(COMET::IReconBase::kTPC1)   ) det=1;
  if (object->UsesDetector(COMET::IReconBase::kTPC2)   ) det=det*10+2;
  if (object->UsesDetector(COMET::IReconBase::kTPC3)   ) det=det*10+3;
  if (object->UsesDetector(COMET::IReconBase::kFGD1)   ) det=det*10+4;
  if (object->UsesDetector(COMET::IReconBase::kFGD2)   ) det=det*10+5;
  if (object->UsesDetector(COMET::IReconBase::kDSECal) ) det=det*10+6;
  if (object->UsesDetector(COMET::IReconBase::kTECal) || 
      object->UsesDetector(COMET::IReconBase::kPECal)  ) det=det*10+7;
  if (object->UsesDetector(COMET::IReconBase::kP0D)    ) det=det*10+8;
  if (object->UsesDetector(COMET::IReconBase::kSMRD)   ) det=det*10+9;

  return det;

}




//*************************************************************
bool COMET::IReconTrackerModule::FillTree(COMET::ICOMETEvent& event) {
//*************************************************************
    //COMETDebug("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!filling the TRACKER tree ");

    fTracks->Clear();

    fNVertices = 0;
    fNTracks = 0;
    fNFGDOther = 0;
    fNTPCOther = 0;
    fNTPCIso = 0;
    fNFGDUnused = 0;
    fNTPCUnused = 0;
    fNTPCExtrapolation = 0;

    fVertices->Clear();
    fTracks->Clear();
    fFGDOther->Clear();
    fTPCOther->Clear();
    fTPCIso->Clear();
    fTPCUnused->Clear();
    fFGDUnused->Clear();
    fTPCExtrapolation->Clear();

    // Get the main result to save.
    COMET::IHandle<COMET::IAlgorithmResult> Result = event.GetFit("ReconTracker");
    if ( !Result ) return true;

    // Now go through each of the different ReconObjects to save:
    COMET::IHandle<COMET::IReconObjectContainer> final = Result->GetResultsContainer("final");         
    COMET::IHandle<COMET::IReconObjectContainer> fgdOther = Result->GetResultsContainer("FGD:other");
    COMET::IHandle<COMET::IReconObjectContainer> tpcOther = Result->GetResultsContainer("TPC:other");
    COMET::IHandle<COMET::IReconObjectContainer> tpcIso = Result->GetResultsContainer("TPC:IsoObjects");
    COMET::IHandle<COMET::IReconObjectContainer> vertices = Result->GetResultsContainer("vertices");
    COMET::IHandle<COMET::IHitSelection> tpcUnused = Result->GetHitSelection("TPC:unused");
    COMET::IHandle<COMET::IHitSelection> fgdUnused = Result->GetHitSelection("FGD:unused");

    // fill fgd other branches
    if ( fgdOther ){
      for ( COMET::IReconObjectContainer::iterator it=fgdOther->begin(); it!=fgdOther->end(); it++ ){
  COMET::IHandle< COMET::IReconTrack > track = (*it);
  if ( track ){
    FillOther( track, fFGDOther, fNFGDOther );
    fNFGDOther++;
  }
      }
    }
    
    // fill tpc other branches
    if ( tpcOther ){
      for ( COMET::IReconObjectContainer::iterator it=tpcOther->begin(); it!=tpcOther->end(); it++ ){
  COMET::IHandle< COMET::IReconTrack > track = (*it);
  if ( track ){
    FillOther( track, fTPCOther, fNTPCOther );
    fNTPCOther++;
  }
      }
    }

    // fill tpc isolated track branches
    // Are there ever any of these?
    if (tpcIso){
      for ( COMET::IReconObjectContainer::iterator it=tpcIso->begin(); it!=tpcIso->end(); it++ ){
  COMET::IHandle< COMET::IReconPID > pid = (*it);
  if ( pid ){
    FillTPC( pid, fTPCIso, fNTPCIso );
    fNTPCIso++;
  }
      }

    }
    
    // fill vertices
    // Again, are there ever any of these?
    if (vertices){
      for ( COMET::IReconObjectContainer::iterator it=vertices->begin(); it!=vertices->end(); it++ ){
  COMET::IHandle< COMET::IReconVertex > vert = (*it);
  if ( vert ){
    FillVertex( vert, fVertices, fNVertices );
    fNVertices++;
  }
      }
    }
    
    // fill unused fgd hits
    if (fgdUnused){
      for ( COMET::IHitSelection::const_iterator hsel = fgdUnused->begin(); hsel != fgdUnused->end(); hsel++){
  IUnusedHit * hit = new ( (*fFGDUnused)[fNFGDUnused++] ) IUnusedHit;
  hit->TotalCharge = (*hsel)->GetCharge();
  hit->Position    = (*hsel)->GetPosition();
  hit->Variance    = (*hsel)->GetUncertainty();
  hit->Time        = (*hsel)->GetTime();
  hit->TimeUnc     = (*hsel)->GetTimeUncertainty();
      }
    }

    // Now loop over the final fitted tracks or pids
    double t0avg = 0.0;
    int NT0 = 0;
    if(final){
      for(COMET::IReconObjectContainer::iterator tt = final->begin(); tt != final->end(); tt++) {
   COMET::IHandle<COMET::IReconPID> pid = *tt;
   COMET::IHandle<COMET::IReconTrack> track = *tt;
   COMET::IHandle<COMET::IReconBase> reco = *tt;
   if ( reco ){
     ITrackerResult * result = new ( (*fTracks)[fNTracks++] ) ITrackerResult;
                 // Save the unique object ID to allow matching to global recon objects
                 result->UniqueID = reco->GetUniqueID();
     result->AlgorithmName = reco->GetAlgorithmName();
     result->Detectors     = GetDetectorNumber( reco );
     result->Status        = reco->CheckStatus( reco->kSuccess );
     result->NDOF          = reco->GetNDOF();
     result->Chi2          = reco->GetQuality();
     if (result->NDOF>0 ) result->Quality = TMath::Prob( result->Chi2 , result->NDOF );
     else result->Quality = -1.0;

     if (pid) {
       COMET::IHandle<COMET::IPIDState> state = pid->GetState();
       if (state){
         t0avg += state->GetPosition().T();
         NT0++;
       }
       FillFinalPid( result, pid );
     }
     else if (track) {
       FillFinalTrack( result, track );
     }

       // Using oaUtility tool to associate the track and the a G4track
       COMET::IG4Trajectory *G4track;
       G4track = NULL;
       double pur = -0xABCDEF;
       double eff = -0xABCDEF;
       int G4TrkId = TrackTruthInfo::GetG4TrajectoryID(reco, eff, pur);
      COMET::IHandle<COMET::IG4TrajectoryContainer> G4trajectories = event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");
      if(G4trajectories){
        COMET::IG4TrajectoryContainer::iterator trajIter = G4trajectories->begin();
        COMET::IG4TrajectoryContainer::iterator trajEnd = G4trajectories->end();
        for( ; trajIter != trajEnd; ++trajIter)
        {
          COMET::IG4Trajectory *trajectory = &(trajIter->second);
          if (trajectory->GetTrackId() == G4TrkId ){
                   G4track = new COMET::IG4Trajectory(*trajectory);
                   break;
          }
        }
             }
             if(G4track!=NULL){
     FillTrueParticle(G4track, pur, eff, result->TrueParticle);
                 }
             }
    }
  }

    if (NT0) t0avg /= double( NT0 );

    // fill unused tpc hits... note we will make use of the average track time
    // to guess the t0 of the unused hits.
    if (tpcUnused){
      COMET::IHandle< COMET::IHitSelection > hstpc = GetHitsFromComboHits( tpcUnused );
      for ( COMET::IHitSelection::const_iterator hsel = hstpc->begin(); hsel != hstpc->end(); hsel++){
  IUnusedHit * hit = new ( (*fTPCUnused)[fNTPCUnused++] ) IUnusedHit;
  //hit->TotalCharge = (*hsel)->GetCharge();
  hit->Position    = (*hsel)->GetPosition();
  hit->Variance    = (*hsel)->GetUncertainty();
  //hit->Time        = (*hsel)->GetTime();
  hit->TimeUnc     = (*hsel)->GetTimeUncertainty();
  
  
  double peakTime=0;
  double maxCharge = 0; 
  
  // Find the pick charge and Time of the Wave form
  COMET::IHandle<COMET::IMultiHit> mh = *hsel;
  std::vector< COMET::IHandle< COMET::ISingleHit > >::const_iterator mhit;
  for( mhit = mh->begin(); mhit !=  mh->end(); mhit++ ) {
    if( maxCharge < (*mhit)->GetCharge() ) {
      maxCharge = (*mhit)->GetCharge(); 
      peakTime = (*mhit)->GetTime();
    }
  }
    
  hit->TotalCharge = maxCharge;

  COMET::IGeometryId geomId = (*hsel)->GetGeomId();
  double t = peakTime - t0avg;
  double x = (fMaxDrift + fCathodeOffset - t*fDriftVelocity)*COMET::IGeomInfo::Get().TPC().GetDriftSense(geomId);

  hit->Position[0] = x;
  hit->Time = t;
      }
    }
 

    // Do the hit association with Claudio's algorithm; this uses the raw TPC result (DR, Sept '10)
    COMET::IHandle<COMET::IAlgorithmResult> TpcResult = event.GetFit("ReconTPC");
    if (!TpcResult)
      return true;

    COMET::IHandle<COMET::IReconObjectContainer> TrackCont = TpcResult->GetResultsContainer("TPCPids");
    COMET::IHandle<COMET::IHitSelection> fgdHits = event.GetHitSelection("fgd");
    if (TrackCont && fgdHits) {
      for (COMET::IReconObjectContainer::iterator it=TrackCont->begin(); it!=TrackCont->end(); it++){
  COMET::IHandle< COMET::IReconPID > pid = *it;
  if (pid) {
    if (ExtrapolateTPC(pid, fTPCExtrapolation, fNTPCExtrapolation, fgdHits))
      fNTPCExtrapolation++;
  }
      }
    }

    return true;
}


COMET::IG4Trajectory COMET::IReconTrackerModule::GetTrajectory(COMET::ICOMETEvent& event,COMET::IHandle<COMET::IReconBase> reco){
  
  COMET::IG4Trajectory traj;
  COMET::IHandle<COMET::IHitSelection> hits = reco->GetHits();

  COMET::IHandle<COMET::IG4TrajectoryContainer> trajectories = event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");
  if (!trajectories) return COMET::IG4Trajectory();
  std::vector<int> search; 

  search.resize(trajectories->size());

  for (COMET::IHitSelection::const_iterator ht = hits->begin(); ht != hits->end(); ht++) { // Loop over reconstructed hits 
    // Get Combo hit
    COMET::IHandle<COMET::IComboHit> ch = *ht;

    // Loop over hits in Combo
    for (COMET::IHitSelection::const_iterator ch_member = ch->GetHits().begin();
   ch_member != ch->GetHits().end(); ++ch_member) {
      COMET::IHandle<COMET::IMultiHit> mh = *ch_member;

      if (!mh) {  // this is the old method with no waveforms.
  std::vector<COMET::IG4VHit*> contributors = HitTruthInfo::GetHitTruthInfo(*ch_member);

  for (std::vector<COMET::IG4VHit*>::const_iterator h = contributors.begin();
       h != contributors.end(); ++h) {
    if (!*h) 
      continue; 
    
    COMET::IG4HitSegment *g4hitseg = dynamic_cast<COMET::IG4HitSegment*>(*h);
    if (!g4hitseg)
      continue;
      
    for (unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
      unsigned int vv = g4hitseg->GetContributors()[i];
      if (vv >= search.size()) 
        search.resize(vv+1);
      search[vv]=search[vv]+1;
    }
  } 
      } else { // This is the new method with waveform
  for (COMET::IMultiHit::iterator wfit = mh->begin(); wfit != mh->end(); wfit++) {
    std::vector<COMET::IG4VHit*> contributors = HitTruthInfo::GetHitTruthInfo(*wfit);	  
    
    for (std::vector<COMET::IG4VHit*>::const_iterator h = contributors.begin();
         h != contributors.end(); ++h) {
      if (!*h) 
        continue;
      
      COMET::IG4HitSegment *g4hitseg = dynamic_cast<COMET::IG4HitSegment*>(*h);
      if (!g4hitseg)
        continue;
      
      for (unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
        unsigned int vv = g4hitseg->GetContributors()[i];
        if (vv >= search.size()) 
      search.resize(vv+1);
        search[vv]=search[vv]+1;
      }
    } 
  }       
      }
    } 
  }

  int maxim = 0; 
  int imax = -1;

  for (unsigned int i = 0; i < search.size() ; i++){
    if (search[i] > maxim) {
      maxim = search[i];
      imax = i; 
    }
  }

  if (imax >= 0) {
    return (*trajectories)[imax];
  } else {
    std::cerr << " TRACK NOT FOUND " << hits->size() << "  " <<  search.size()  <<  std::endl; 
    return COMET::IG4Trajectory();
  }
}



//*************************************************************
void COMET::IReconTrackerModule::FillTrueParticle( COMET::IG4Trajectory* G4track, double pur, double eff, COMET::IReconTrackerModule::ITrueParticle& trueParticle ) {
//*************************************************************

  if (G4track){
    
    trueParticle.ID    = G4track->GetTrackId();
    trueParticle.Pur = pur;
    trueParticle.Eff = eff;

    COMET::IG4PrimaryVertex G4vertex;
    bool ok = GetG4Vertex(G4track,G4vertex);
    FillTrueVertex(ok, G4vertex, 0, 0, trueParticle.Vertex);

  } else {
    
    //std::cout<<"True G4 track not found"<<std::endl;
    trueParticle.ID    = -1;

    COMET::IG4PrimaryVertex G4vertex;
    FillTrueVertex(false, G4vertex, 0, 0, trueParticle.Vertex);
  }
  return;
}

//*************************************************************
COMET::IHandle<COMET::IG4Trajectory> COMET::IReconTrackerModule::GetParent(COMET::IG4Trajectory* G4track) {
//*************************************************************

  const COMET::ICOMETEvent& event = *(COMET::IEventFolder::GetCurrentEvent());  
  COMET::IHandle<COMET::IG4TrajectoryContainer> trajectories = event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");

  if (G4track->GetParentId()!=0){ 
    COMET::IHandle<COMET::IG4Trajectory> traj(new COMET::IG4Trajectory( ( (*trajectories)[G4track->GetParentId()]  )));
    return traj;
  }				
  else
    return COMET::IHandle<COMET::IG4Trajectory>();

}

//*************************************************************
COMET::IHandle<COMET::IG4Trajectory> COMET::IReconTrackerModule::GetParent(COMET::IHandle<IG4Trajectory> G4track) {
//*************************************************************

  const COMET::ICOMETEvent& event = *(COMET::IEventFolder::GetCurrentEvent());  
  COMET::IHandle<COMET::IG4TrajectoryContainer> trajectories = event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");

  if (G4track->GetParentId()!=0){ 
    COMET::IHandle<COMET::IG4Trajectory> traj(new COMET::IG4Trajectory( ( (*trajectories)[G4track->GetParentId()]  )));
    return traj;
  }				
  else
    return COMET::IHandle<COMET::IG4Trajectory>();

}

//*************************************************************
void COMET::IReconTrackerModule::FillTrueVertex( bool found, const COMET::IG4PrimaryVertex& G4vertex, double pur, double eff, COMET::ITrueVertex& vertex) {
//*************************************************************

  if (found){
    vertex.Position = G4vertex.GetPosition();
    vertex.ID = G4vertex.GetInteractionNumber();
  }
  else{
    vertex.Position = TLorentzVector(0,0,0,0);
    vertex.ID = -1;
  }

}
  

//*************************************************************
bool  COMET::IReconTrackerModule::GetG4Vertex(COMET::IG4Trajectory* G4track, COMET::IG4PrimaryVertex&  G4vertex){
//*************************************************************

  const COMET::ICOMETEvent& event = *(COMET::IEventFolder::GetCurrentEvent());
  COMET::IHandle<COMET::IG4PrimaryVertexContainer>   vertices  = event.Get<COMET::IG4PrimaryVertexContainer>("truth/G4PrimVertex00");

  if (vertices){
    std::vector<COMET::IG4PrimaryVertex>::iterator it1;
    for(it1= vertices->begin(); it1!= vertices->end(); it1++){
      const COMET::TG4PrimaryParticleContainer& particles = it1->GetPrimaryParticles();
      
      std::vector<COMET::IG4PrimaryParticle>::const_iterator it2;
      for(it2 = particles.begin(); it2 != particles.end(); it2++) {
  if (it2->GetTrackId() == G4track->GetTrackId()){
    G4vertex = *it1;
    return true;
  }
      }
    }
  }

  return false;

}


//*************************************************************
bool  COMET::IReconTrackerModule::GetIncomingParticle(const COMET::IG4PrimaryVertex&  G4vertex, COMET::IG4PrimaryParticle& incoming){
//*************************************************************

  for(std::vector<COMET::IG4PrimaryVertex>::const_iterator infoVtxIter = G4vertex.GetInfoVertex().begin();
      infoVtxIter != G4vertex.GetInfoVertex().end(); ++infoVtxIter) {
    std::vector<COMET::IG4PrimaryParticle>::const_iterator particleIter;
    for(particleIter = infoVtxIter->GetPrimaryParticles().begin();
  particleIter != infoVtxIter->GetPrimaryParticles().end(); ++particleIter) {
      Int_t pdg = particleIter->GetPDGCode();
      if(particleIter->GetTrackId() == -1) {
  // is an incoming particle (including target)
  if(TMath::Abs(pdg) == 12 || TMath::Abs(pdg) == 14 || TMath::Abs(pdg) == 16) {
    incoming =  *particleIter;
    return true;
  }
      }    
    }
  }

  return false;

}


//*************************************************************
int COMET::IReconTrackerModule::GetNumberOfHits(COMET::IHandle<COMET::IReconBase> object) {
//*************************************************************

  int nHits = 0;
  if (object->GetHits())
    nHits += object->GetHits()->size();
  else if (object->GetConstituents()){
    COMET::IReconObjectContainer::const_iterator it;
    for (it=object->GetConstituents()->begin(); it != object->GetConstituents()->end(); it++){
      COMET::IHandle<COMET::IReconBase> obj = (*it);
      if (obj ){
  nHits += GetNumberOfHits( obj );
      }
    }
  }

  return nHits;

}