#include "IGlobalReconModule.hxx"

#include <functional>
#include <algorithm>

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

#include <IGeomInfo.hxx>

#include <IComboHit.hxx>
#include <IReconHit.hxx>
#include <IG4HitSegment.hxx>
#include <ITPCmcDigit.hxx>

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

#include "fgdUtils.hxx"
#include "IIntegerDatum.hxx"
#include <TrackTruthInfo.hxx>
#include <TrackingUtils.hxx>
#include <ReconObjectUtils.hxx>

const int NMAXALTERNATES = 3;
const int NMAXECAL = 3;
const int NMAXP0DECAL = 2;
const int NMAXTPC  = 3;
const int NMAXFGD  = 2;
const int NMAXP0D  = 2;
const int NMAXSMRD = 3;
const int NMAXTRACKER = 10;

const int NMAXTPCOTHER  = 20;

const bool debug_kine = false;
const bool debug_pids = false;
const bool debug_broken = false;
const bool debug_extrap = false;
const bool debug_vertex_extrap = false;

const int DEFAULT_MIN =  99999;
const int DEFAULT_MAX = -99999;

ClassImp(COMET::IGlobalReconModule);

ClassImp(COMET::IGlobalReconModule::ITpcPID);
COMET::IGlobalReconModule::ITpcPID::~ITpcPID(){}

ClassImp(COMET::IGlobalReconModule::IGlobalPID);
COMET::IGlobalReconModule::IGlobalPID::~IGlobalPID(){/*if (TrueParticle) delete TrueParticle;*/}

COMET::IGlobalReconModule::IGlobalPID::IGlobalPID(){
  Alternates = new TClonesArray("COMET::IGlobalReconModule::IGlobalPIDAlternate",  NMAXALTERNATES );
  ECAL = new TClonesArray("COMET::IGlobalReconModule::IECALObject", NMAXECAL);
  P0DECAL = new TClonesArray("COMET::IGlobalReconModule::IECALObject", NMAXP0DECAL);
  TPC  = new TClonesArray("COMET::IGlobalReconModule::ITPCObject",  NMAXTPC );
  FGD  = new TClonesArray("COMET::IGlobalReconModule::IFGDObject",  NMAXFGD );
  P0D  = new TClonesArray("COMET::IGlobalReconModule::IP0DObject",  NMAXP0D );
  SMRD = new TClonesArray("COMET::IGlobalReconModule::ISMRDObject", NMAXSMRD);
  TRACKER  = new TClonesArray("COMET::IGlobalReconModule::ITrackerObject",  NMAXTRACKER );
  HitsSaved  = new TClonesArray("COMET::IGlobalReconModule::IGlobalHit",  4 );

  // TPC digit decoding is a little too verbose. 
  // Disable it here.
  COMET::ICOMETLog::SetLogLevel("oaUnpack", COMET::ICOMETLog::QuietLevel);


  NAlternates=0;
  Alternates->Clear();

  NDsECALs=0;
  NTrECALs=0;
  NECALs=0;
  ECAL->Clear();


  NP0DECALs=0;
  P0DECAL->Clear();


  NTPCs=0;
  TPC->Clear();

  NFGDs=0;
  FGD->Clear();

  NSMRDs=0;
  SMRD->Clear();

  NP0Ds=0;
  P0D->Clear();

  NTRACKERs=0;
  TRACKER->Clear();

  NHitsSaved=0;
  HitsSaved->Clear();


}

ClassImp(COMET::IGlobalReconModule::IVertexConstituent);
COMET::IGlobalReconModule::IVertexConstituent::~IVertexConstituent(){}

ClassImp(COMET::IGlobalReconModule::IGlobalVertex);
COMET::IGlobalReconModule::IGlobalVertex::~IGlobalVertex(){}
COMET::IGlobalReconModule::IGlobalVertex::IGlobalVertex(){

  Constituents  = new TClonesArray("COMET::IGlobalReconModule::IVertexConstituent",  NCONSTITUENTS);
  TrueVertices = new TClonesArray("COMET::ITrueVertex", NCONSTITUENTS);
  NConstituents = 0;
  Constituents->Clear();
  NTrueVertices = 0;
  TrueVertices->Clear();
}

//ClassImp(COMET::IGlobalReconModule::ISubBaseObject);
//COMET::IGlobalReconModule::ISubBaseObject::~ISubBaseObject(){}

ClassImp(COMET::IGlobalReconModule::IP0DObject);
COMET::IGlobalReconModule::IP0DObject::~IP0DObject(){}

ClassImp(COMET::IGlobalReconModule::IECALObject);
COMET::IGlobalReconModule::IECALObject::~IECALObject(){}

ClassImp(COMET::IGlobalReconModule::ITPCObject);
COMET::IGlobalReconModule::ITPCObject::~ITPCObject(){}

ClassImp(COMET::IGlobalReconModule::ITPCOtherObject);
COMET::IGlobalReconModule::ITPCOtherObject::~ITPCOtherObject(){}


ClassImp(COMET::IGlobalReconModule::IFGDObject);
COMET::IGlobalReconModule::IFGDObject::~IFGDObject(){}

ClassImp(COMET::IGlobalReconModule::ISMRDObject);
COMET::IGlobalReconModule::ISMRDObject::~ISMRDObject(){}

ClassImp(COMET::IGlobalReconModule::ITrackerObject);
COMET::IGlobalReconModule::ITrackerObject::~ITrackerObject(){}

ClassImp(COMET::IGlobalReconModule::IFgdCluster);
COMET::IGlobalReconModule::IFgdCluster::~IFgdCluster(){}

ClassImp(COMET::IGlobalReconModule::IFgdTimeBin);
COMET::IGlobalReconModule::IFgdTimeBin::~IFgdTimeBin(){}

COMET::IGlobalReconModule::IFgdTimeBin::IFgdTimeBin(){

  FGD1Unused = new TClonesArray("COMET::IGlobalReconModule::IGlobalHit", 20);
  FGD2Unused = new TClonesArray("COMET::IGlobalReconModule::IGlobalHit", 20);

  NFGD1Unused=0;
  FGD1Unused->Clear();

  NFGD2Unused=0;
  FGD2Unused->Clear();

}


ClassImp(COMET::IGlobalReconModule::IGlobalHit);
COMET::IGlobalReconModule::IGlobalHit::~IGlobalHit(){}

ClassImp(COMET::IGlobalReconModule::ISMRDHit);
COMET::IGlobalReconModule::ISMRDHit::~ISMRDHit(){}


ClassImp(COMET::IGlobalReconModule::IOutermostHits);
COMET::IGlobalReconModule::IOutermostHits::~IOutermostHits(){}


bool SortHitsInX(const COMET::IHandle<COMET::IHit>& h1, const COMET::IHandle<COMET::IHit>& h2){
  return h1->GetPosition().X() < h2->GetPosition().X();
}
bool SortHitsInY(const COMET::IHandle<COMET::IHit>& h1, const COMET::IHandle<COMET::IHit>& h2){
  return h1->GetPosition().Y() < h2->GetPosition().Y();
}
bool SortHitsInZ(const COMET::IHandle<COMET::IHit>& h1, const COMET::IHandle<COMET::IHit>& h2){
  return h1->GetPosition().Z() < h2->GetPosition().Z();
}


bool SortNodesInX(const COMET::IHandle<COMET::IReconNode>& n1, const COMET::IHandle<COMET::IReconNode>& n2){
  return TrackingUtils::GetPosition(n1->GetState()).X() < TrackingUtils::GetPosition(n2->GetState()).X();
}
bool SortNodesInY(const COMET::IHandle<COMET::IReconNode>& n1, const COMET::IHandle<COMET::IReconNode>& n2){
  return TrackingUtils::GetPosition(n1->GetState()).Y() < TrackingUtils::GetPosition(n2->GetState()).Y();
}
bool SortNodesInZ(const COMET::IHandle<COMET::IReconNode>& n1, const COMET::IHandle<COMET::IReconNode>& n2){
  return TrackingUtils::GetPosition(n1->GetState()).Z() < TrackingUtils::GetPosition(n2->GetState()).Z();
}
bool SortNodesInXReverse(const COMET::IHandle<COMET::IReconNode>& n1, const COMET::IHandle<COMET::IReconNode>& n2){
  return TrackingUtils::GetPosition(n1->GetState()).X() >= TrackingUtils::GetPosition(n2->GetState()).X();
}
bool SortNodesInYReverse(const COMET::IHandle<COMET::IReconNode>& n1, const COMET::IHandle<COMET::IReconNode>& n2){
  return TrackingUtils::GetPosition(n1->GetState()).Y() >= TrackingUtils::GetPosition(n2->GetState()).Y();
}
bool SortNodesInZReverse(const COMET::IHandle<COMET::IReconNode>& n1, const COMET::IHandle<COMET::IReconNode>& n2){
  return TrackingUtils::GetPosition(n1->GetState()).Z() >= TrackingUtils::GetPosition(n2->GetState()).Z();
}


#define CVSTAG "\
  $Name: v5r19 $"
#define CVSID  "\
  $Id: IGlobalReconModule.cxx,v 1.122 2013/02/22 15:03:27 db211 Exp $"

//*************************************************************
COMET::IGlobalReconModule::IGlobalReconModule(const char *name, const char *title){
//*************************************************************

  fTestTPCInfo = false;

  SetNameTitle(name, title);

  // Enable this module by default:
  fIsEnabled = kTRUE;
  fDescription = "Global Recon Output";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;

  // Tree variables
  fNPIDs = 0;
  fPIDs = new TClonesArray("COMET::IGlobalReconModule::IGlobalPID", 20);

  fNTPCOthers = 0;
  fTPCOthers = new TClonesArray("COMET::IGlobalReconModule::ITPCOtherObject", NMAXTPCOTHER);

  fNVertices = 0;
  fVertices = new TClonesArray("COMET::IGlobalReconModule::IGlobalVertex", 20);
  fPVInd = 0;
  fNSVertices = 0;

  fNDelayedClusters = 0;
  fDelayedClusters = new TClonesArray("COMET::IGlobalReconModule::IFgdCluster", 10);

  fNFgdTimeBins = 0;
  fFgdTimeBins = new TClonesArray("COMET::IGlobalReconModule::IFgdTimeBin", 20);

  fSMRDUnused = new TClonesArray("COMET::IGlobalReconModule::ISMRDHit", 100);

  fP0DUnused = new TClonesArray("COMET::IGlobalReconModule::IGlobalHit", 100);

  if (fTestTPCInfo){
    fNTPCPIDs = 0;
    fTPCPIDs = new TClonesArray("COMET::IGlobalReconModule::ITpcPID", 20);
  }
  
  fRecPackInitialized = false;



}

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

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

  return true;
}

//*************************************************************
void COMET::IGlobalReconModule::InitializeModule() {}
//*************************************************************

//*************************************************************
void COMET::IGlobalReconModule::InitializeBranches() {
//*************************************************************
  SetSplitLevel(1);


  if (fTestTPCInfo){
    fOutputTree->Branch("NTPCPIDs",
      &fNTPCPIDs, "NTPCPIDs/I",
      fBufferSize);
    
    fOutputTree->Branch("TPCPIDs","TClonesArray",
      &fTPCPIDs,
      fBufferSize, fSplitLevel);
  }


  // The number of vertex objects that have been saved.
  fOutputTree->Branch("NVertices",
          &fNVertices, "NVertices/I",
          fBufferSize);

  //Vertex branch
  fOutputTree->Branch("Vertices",
          &fVertices,
          fBufferSize, fSplitLevel);

  // Delayed clusters in FGD
  fOutputTree->Branch("NDelayedClusters",
          &fNDelayedClusters, "NDelayedClusters/I",
          fBufferSize);


  fOutputTree->Branch("DelayedClusters",
          &fDelayedClusters,
          fBufferSize, fSplitLevel);
 

  // FGD hit time bins
  fOutputTree->Branch("NFgdTimeBins",
          &fNFgdTimeBins, "NFgdTimeBins/I",
          fBufferSize);


  fOutputTree->Branch("FgdTimeBins",
          &fFgdTimeBins,
          fBufferSize, fSplitLevel);
 

  // The number of objects that have been saved.
  fOutputTree->Branch("NPIDs",
          &fNPIDs, "NPIDs/I",
          fBufferSize);

  // A branch for the clone array of Global object information.
  fOutputTree->Branch("PIDs","TClonesArray",
          &fPIDs,
          fBufferSize, fSplitLevel);


  // The number of objects that have been saved.
  fOutputTree->Branch("NTPCOthers",
          &fNTPCOthers, "NTPCOthers/I",
          fBufferSize);

  // A branch for the clone array of Global object information.
  fOutputTree->Branch("TPCOthers","TClonesArray",
          &fTPCOthers,
          fBufferSize, fSplitLevel);


  //number of unused hits in FGD1, FGD2, P0D and SMRD
  fOutputTree->Branch("NFGD1Unused", &fNFGD1Unused, "NFGD1Unused/I", fBufferSize);                       
  fOutputTree->Branch("NFGD2Unused", &fNFGD2Unused, "NFGD2Unused/I", fBufferSize);                       


  fOutputTree->Branch("NP0DUnused", &fNP0DUnused, "NP0DUnused/I", fBufferSize);                       


  fOutputTree->Branch("P0DUnused","TClonesArray",
          &fP0DUnused,
          fBufferSize, fSplitLevel);
             
  fOutputTree->Branch("SMRDUnused","TClonesArray",&fSMRDUnused,fBufferSize, fSplitLevel); 
  fOutputTree->Branch("NSMRDUnused", &fNSMRDUnused, "NSMRDUnused/I", fBufferSize);  
  fOutputTree->Branch("NSMRDTopUnused",    &fNSMRDTopUnused,    "NSMRDTopUnused/I", fBufferSize);                       
  fOutputTree->Branch("NSMRDBottomUnused", &fNSMRDBottomUnused, "NSMRDBottomUnused/I", fBufferSize);                       
  fOutputTree->Branch("NSMRDLeftUnused",   &fNSMRDLeftUnused,   "NSMRDLeftUnused/I", fBufferSize);                       
  fOutputTree->Branch("NSMRDRightUnused",  &fNSMRDRightUnused,  "NSMRDRightUnused/I", fBufferSize);                       


  fOutputTree->Branch("NTPCUnused", &fNTPCUnused, "NTPCUnused/I", fBufferSize);                       

  //Save P0D outermost hits
  fOutputTree->Branch("P0DOutermostHits", "IOutermostHits",&fP0DOutermostHits, fBufferSize,fSplitLevel);                       


  //time of the median hit of the earliest track
  fOutputTree->Branch("EarliestTrackTime", &fEarliestTrackMedianHitTime, fBufferSize);  


  // Counters for each detector
  fOutputTree->Branch("NDsECAL",  &fNDsECAL,  "NDsECAL/I", fBufferSize);
  fOutputTree->Branch("NTrECAL",  &fNTrECAL,  "NTrECAL/I", fBufferSize);
  fOutputTree->Branch("NP0DECAL", &fNP0DECAL, "NP0DECAL/I", fBufferSize);
  fOutputTree->Branch("NTPC",  &fNTPC,  "NTPC/I",  fBufferSize);
  fOutputTree->Branch("NFGD",  &fNFGD,  "NFGD/I",  fBufferSize);
  fOutputTree->Branch("NP0D",  &fNP0D,  "NP0D/I",  fBufferSize);
  fOutputTree->Branch("NSMRD", &fNSMRD, "NSMRD/I", fBufferSize);

  fOutputTree->Branch("NDsECALIso",  &fNDsECALIso,  "NDsECALIso/I", fBufferSize);
  fOutputTree->Branch("NTrECALIso",  &fNTrECALIso,  "NTrECALIso/I", fBufferSize);
  fOutputTree->Branch("NP0DECALIso", &fNP0DECALIso, "NP0DECALIso/I", fBufferSize);
  fOutputTree->Branch("NTPCIso",  &fNTPCIso,  "NTPCIso/I",  fBufferSize);
  fOutputTree->Branch("NFGDIso",  &fNFGDIso,  "NFGDIso/I",  fBufferSize);
  fOutputTree->Branch("NP0DIso",  &fNP0DIso,  "NP0DIso/I",  fBufferSize);
  fOutputTree->Branch("NSMRDIso", &fNSMRDIso, "NSMRDIso/I", fBufferSize);


}



//*************************************************************
bool COMET::IGlobalReconModule::FillTree(COMET::ICOMETEvent& event) {
//*************************************************************
    
  if(!fRecPackInitialized){
    // add a new RecpackManager 
    COMET::tman().add_manager("globalAnalysis");
    COMET::tman().select_manager("globalAnalysis");
    
    // only surfaces with measurements are tried inside a volume
    COMET::rpman("globalAnalysis").navigation_svc().navigator(RP::particle_helix).set_unique_surface(true);
    
    // enable energy loss correction
    COMET::rpman("globalAnalysis").model_svc().enable_correction(RP::particle_helix, RP::eloss, true);
    
    //enable MS
    COMET::rpman("globalAnalysis").model_svc().enable_noiser(RP::particle_helix, RP::ms, true);
    
    fEquation = COMET::rpman().model_svc().model(RP::particle_helix).retrieve_tool<HelixEquation>("equation");

    if (debug_extrap){
      COMET::tman("globalAnalysis").SetVerbosity(0,2,0,0);
      COMET::rpman("globalAnalysis").model_svc().model(RP::particle_helix).tool("correction/eloss").set_verbosity(Messenger::VVERBOSE);
      COMET::rpman("globalAnalysis").navigation_svc().navigator(RP::particle_helix).master_inspector().set_verbosity(Messenger::VVERBOSE);
      COMET::rpman("globalAnalysis").navigation_svc().inspector("eloss").set_verbosity(Messenger::VVERBOSE);     
    }
    InitializeExtrapolationToDetectors();

    fRecPackInitialized = true;
  }
  
  if (fTestTPCInfo){
    fNTPCPIDs = 0;
    fTPCPIDs->Clear();
  }

  fNVertices = 0;
  fVertices->Clear();
  fNPIDs = 0;
  fNSVertices = 0;

  fNPIDs = 0;
  fPIDs->Clear();

  fNTPCOthers = 0;
  fTPCOthers->Clear();

  fNDelayedClusters = 0;
  fDelayedClusters->Clear();

  fNFgdTimeBins = 0;
  fFgdTimeBins->Clear();

  fNFGD1Unused = 0;
  fNFGD2Unused = 0;

  fNP0DUnused = 0;
  fP0DUnused->Clear();

  fSMRDUnused->Clear(); 
  fNSMRDUnused = 0;
  fNSMRDTopUnused = 0;
  fNSMRDBottomUnused = 0;
  fNSMRDLeftUnused = 0;
  fNSMRDRightUnused = 0;

  fNTPCUnused = 0;

  fEarliestTrackMedianHitTime = DEFAULT_MIN; 

  fNDsECAL = 0;
  fNTrECAL = 0;
  fNP0DECAL = 0;
  fNTPC = 0;
  fNFGD = 0;
  fNP0D = 0;
  fNSMRD = 0;

  fNDsECALIso = 0;
  fNTrECALIso = 0;
  fNP0DECALIso = 0;
  fNTPCIso = 0;
  fNFGDIso = 0;
  fNP0DIso = 0;
  fNSMRDIso = 0;


  //  std::cout << "event ID: " << event.GetEventId() << std::endl;

  // Get the main result to save.
  COMET::IHandle<COMET::IAlgorithmResult> globalRecon = event.GetFit("ReconGlobal");
  if (!globalRecon){
    //    std::cout << "no ReconGlobal result" << std::endl;
    return true;
  }

  // Get the results from reconstruction
  COMET::IHandle<COMET::IReconObjectContainer> allObjects
    = globalRecon->GetResultsContainer("final");

  if (!allObjects) return true;

  // This is only for tests
  if (fTestTPCInfo){
    COMET::IHandle<COMET::IAlgorithmResult> ReconTPC = event.GetFit("ReconTPC");
    if (ReconTPC){
      
      // Get the results from reconstruction
      COMET::IHandle<COMET::IReconObjectContainer> tpcObjects
  = ReconTPC->GetResultsContainer("TPCPids");
      
      if (tpcObjects){
  // Now loop over objects.
  for(COMET::IReconObjectContainer::iterator tt = tpcObjects->begin();
      tt != tpcObjects->end(); tt++) {
    
    COMET::IHandle<COMET::IReconBase> object = *tt;
    if(object) FillTPCPID(object);
  }
      }
    }
  }

  // fill the map of broken tracks (to be used later)
  FillBrokenTracksMap(*allObjects);

  // Now loop over object and fill the global analysis objects
  fGlobalIndexMap.clear(); 
  for(COMET::IReconObjectContainer::iterator tt = allObjects->begin();
      tt != allObjects->end(); tt++) {
    
    COMET::IHandle<COMET::IReconBase> object = *tt;
    if (!object) continue;

    //The way to be done in future
    //   COMET::IHandle<COMET::IReconVertex> vertex = *tt;      
    //     if (vertex)
    //         FillVertex(object);
    //     else//Fill track/pid informtaion

    FillGlobalPIDs(event, object);
    
  }



  FillVertexInfo(event, allObjects);
 
  FillDelayedFgdClusters(event, allObjects);

  FillFgdTimeBins(event);

  FillUnusedHits(event);

  FillTPCOther(event);

  // Fill P0D outermost hits
  double p0d_charge_cut = 7;  
  COMET::IHandle<COMET::IHitSelection> p0dHits = event.GetHitSelection("p0d");
  if (p0dHits)
    FillOutermostHits(*p0dHits, p0d_charge_cut, fP0DOutermostHits);

  fGlobalIndexMap.clear(); 
  fBrokenIndexMap.clear();

  return true;
}

//*************************************************************
void COMET::IGlobalReconModule::FillGlobalPIDs(COMET::ICOMETEvent& event, COMET::IHandle<COMET::IReconBase> object) {
//*************************************************************

  if (debug_pids)
    std::cout << "Fill Global PID: " << object << std::endl;

  // Get the compact detector number
  int dets = GetDetectorNumber(object);

  // update counters for specfic detectors
  if (object->UsesDetector(COMET::IReconBase::kDSECal) ){
    fNDsECAL++;
    if (dets < 10) fNDsECALIso++;
  }
  if (object->UsesDetector(COMET::IReconBase::kECal)){
    fNTrECAL++;
    if (dets < 10) fNTrECALIso++;
  }
  if (object->UsesDetector(COMET::IReconBase::kPECal)){
    fNP0DECAL++;
    if (dets < 10) fNP0DECALIso++;
  }
  if (object->UsesDetector(COMET::IReconBase::kFGD)){
    fNFGD++;
    if (dets < 10) fNFGDIso++;
  }
  if (object->UsesDetector(COMET::IReconBase::kTPC)){
    fNTPC++;
    if (dets < 10) fNTPCIso++;
  }
  if (object->UsesDetector(COMET::IReconBase::kP0D)){
    fNP0D++;
    if (dets < 10) fNP0DIso++;
  }
  if (object->UsesDetector(COMET::IReconBase::kSMRD)){
    fNSMRD++;
    if (dets < 10) fNSMRDIso++;
  }


  // __________________________________________________________
  // Fill the reconstructed Global PIDs.

  // try to get the state of the object. If it dosen't exists skip the object. 
  try{COMET::IHandle<COMET::IReconState> tstate = object->GetState(); }
  catch(COMET::EReconObject e) { std::cout << "IReconBase with no state. Skip object !!" << std::endl; return; }

  // Fill the map of indexes used for vertex link
  fGlobalIndexMap[object] = fNPIDs;

  // create a new IGlobalPID
  IGlobalPID* globalObject = new((*fPIDs)[fNPIDs++]) IGlobalPID;

  // the unique ID of the global object
  globalObject->UniqueID = object->GetUniqueID();

  //--------- Fill the vector of unique IDs of associated tracks
  if (fBrokenIndexMap.find(object) != fBrokenIndexMap.end())
    globalObject->BrokenUniqueIDs = fBrokenIndexMap[object];

  //--------- Information about the detectors that used in this object
  globalObject->Detectors  = dets;
  FillDetectorUsed(object, globalObject->DetectorUsed);

  //------- General information about the track --------
  
  globalObject->isForward     = (TrackingUtils::GetDirection(object).Z() > 0);  
  globalObject->SenseOK       = TrackingUtils::GetSenseOK(*object);  
  globalObject->Charge        = TrackingUtils::GetCharge(object);	
  globalObject->EDeposit      = TrackingUtils::GetEDeposit(object);
  globalObject->Width         = TrackingUtils::GetWidth(object).X(); //! take only the first component
  globalObject->Cone          = TrackingUtils::GetCone(object);
  globalObject->Length        = ComputeTrackLength(object);  
      
  globalObject->AlgorithmName = object->GetAlgorithmName();
  globalObject->Status     = object->CheckStatus(object->kSuccess);
  globalObject->Chi2       = object->GetQuality();
  globalObject->NDOF       = object->GetNDOF();
  try{ globalObject->NNodes     = object->GetNodes().size();}
  catch(COMET::EReconObject e) { globalObject->NNodes = 0; }
  
  globalObject->NHits = DEFAULT_MAX;  
  if (object->GetHits())
    globalObject->NHits    = object->GetHits()->size();

  // The number of constituents. For composite objects
  globalObject->NConstituents = 0;
  if (object->GetConstituents())
    globalObject->NConstituents  = object->GetConstituents()->size();


  //---------- Fill info about first and last hits
  if (object->GetHits()){
    FillFirstLastHits(*(object->GetHits()), *globalObject);
  }

  //--------- Fill the true particle when it is found
  double pur, eff;
  COMET::IHandle<COMET::IG4Trajectory> G4track = GetG4Trajectory(*object,pur,eff);
  FillTrueParticle(G4track, pur, eff, globalObject->TrueParticle);

  //--------- Fill the kinematic extrapolated to the true vertex
  FillKinematicsAtTrueVertex(G4track, object,
           globalObject->PositionAtTrueVertex,  globalObject->PositionVarAtTrueVertex,
           globalObject->DirectionAtTrueVertex, globalObject->DirectionVarAtTrueVertex,
           globalObject->MomentumAtTrueVertex,  globalObject->MomentumErrorAtTrueVertex);

  //------- Front and Back kinematics    
  COMET::IHandle<COMET::IReconState> firstState = TrackingUtils::GetFirstState(*object);
  FillKinematics(firstState,
     globalObject->FrontPosition, globalObject->FrontPositionVar,
     globalObject->FrontDirection,globalObject->FrontDirectionVar,
     globalObject->FrontMomentum, globalObject->FrontMomentumError);
  
  COMET::IHandle<COMET::IReconState> lastState  = TrackingUtils::GetLastState(*object);
  FillKinematics(lastState,
     globalObject->BackPosition, globalObject->BackPositionVar,
     globalObject->BackDirection,globalObject->BackDirectionVar,
     globalObject->BackMomentum, globalObject->BackMomentumError);

  //------- Momentumby range

  //globalObject->RangeMomentum  = TrackingUtils::GetMomentumByRange(*object);
  globalObject->RangeMomentumMuon = TrackingUtils::GetTRealDatumValue(*object,"prange_muon",DEFAULT_MAX);
  globalObject->RangeMomentumMuonFlip = TrackingUtils::GetTRealDatumValue(*object,"prange_muon_flipped",DEFAULT_MAX);
  globalObject->RangeMomentumElectron = TrackingUtils::GetTRealDatumValue(*object,"prange_electron",DEFAULT_MAX);
  globalObject->RangeMomentumElectronFlip = TrackingUtils::GetTRealDatumValue(*object,"prange_electron_flipped",DEFAULT_MAX);
  globalObject->RangeMomentumProton = TrackingUtils::GetTRealDatumValue(*object,"prange_proton",DEFAULT_MAX);
  globalObject->RangeMomentumProtonFlip = TrackingUtils::GetTRealDatumValue(*object,"prange_proton_flipped",DEFAULT_MAX);



  //----------  Set the PID information
  COMET::IHandle<COMET::IReconPID> PID = object;
  if (PID){      
    globalObject->PIDWeights.push_back(PID->GetPIDWeight());
    globalObject->ParticleIds.push_back( ComputeParticleId(PID) );

    COMET::IReconObjectContainer::const_iterator it;
    for (it=PID->GetAlternates().begin(); it!=PID->GetAlternates().end(); it++){
      COMET::IHandle<COMET::IReconPID> alter = *it;
      globalObject->PIDWeights.push_back(alter->GetPIDWeight());
      globalObject->ParticleIds.push_back( ComputeParticleId(alter) );
    }
  }

  //---------- Extrapolate to the entrance and exit of each subdetector and save position, direction and momentum
  FillExtrapolationToDetectors(object, *globalObject);
  
  //---------- get FGD vertex activity before calling FillFGDInfo but after front/back position defined
  // Evaluate vertex activity at track front position if in either FGD, otherwise use back position
  if( COMET::IGeomInfo::Get().FGD().IsInFGD1( globalObject->FrontPosition.Vect() ) || COMET::IGeomInfo::Get().FGD().IsInFGD2( globalObject->FrontPosition.Vect() ) )
    GetFGDSimpleVA(event, object, globalObject->FrontPosition);
  else if( COMET::IGeomInfo::Get().FGD().IsInFGD1( globalObject->BackPosition.Vect() ) || COMET::IGeomInfo::Get().FGD().IsInFGD2( globalObject->BackPosition.Vect() ) )
    GetFGDSimpleVA(event, object, globalObject->BackPosition);

  //---------- Fill subdetector information
  FillECALInfo(object, *globalObject);
  FillP0DInfo( object, *globalObject);
  FillTPCInfo( object, *globalObject);
  FillFGDInfo( object, *globalObject);
  FillSMRDInfo(object, *globalObject);
  FillTrackerInfo( object, *globalObject);

  //---------- Fill the alternates information
  FillGlobalPIDAlternates(G4track, object, *globalObject);


  // Add subdetector recomputation of some variables for isolated objects
  // For ECAL
  if (globalObject->Detectors<10 && globalObject->NECALs==1){
    IECALObject* sub = (IECALObject*)(*(globalObject->ECAL))[0];
    globalObject->Cone     = sub->Cone;
    globalObject->EDeposit = sub->EDeposit;
    globalObject->Width    = sub->Width.X();
    if (globalObject->Length==0)
      globalObject->Length = sub->Length;
  }
  // For P0D
  else if (globalObject->Detectors<10 && globalObject->NP0Ds==1){
    IP0DObject* sub = (IP0DObject*)(*(globalObject->P0D))[0];
    globalObject->Cone      = sub->Cone;
    globalObject->EDeposit  = sub->EDeposit;
    globalObject->Width     = sub->Width;
    if (globalObject->Length==0)
      globalObject->Length = sub->Length;
  }
}

//*************************************************************
int COMET::IGlobalReconModule::ComputeParticleId(COMET::IHandle<COMET::IReconPID> PID){
//*************************************************************

  // only three hypothesis considered at recon level 

  int pdg=0;

  if      (PID->GetParticleId() == COMET::IReconPID::kElectron) pdg = 11;
  else if (PID->GetParticleId() == COMET::IReconPID::kMuon)     pdg = 13;
  else if (PID->GetParticleId() == COMET::IReconPID::kProton)   pdg = 2212;
    
  return pdg;
}

//*************************************************************
void COMET::IGlobalReconModule::FillKinematics(COMET::IHandle<COMET::IReconState> state,
              TLorentzVector& pos, TLorentzVector& posVar, 
              TVector3& dir, TVector3& dirVar, 
              double& mom, double& momErr){
//*************************************************************

  if (state){
    pos    = TrackingUtils::GetPosition(state);
    posVar = TrackingUtils::GetPositionVariance(state);     
    dir    = TrackingUtils::GetDirection(state);     
    dirVar = TrackingUtils::GetDirectionVariance(state);     
    mom    = TrackingUtils::GetMomentum(state);
    momErr = TrackingUtils::GetMomentumError(state);
  }
  else{
    pos    = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    posVar = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    dir    = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    dirVar = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    mom    = DEFAULT_MAX; 
    momErr  = DEFAULT_MAX; 
  }

  if (debug_kine){
    std::cout << "position:  " << pos.X() << " " << pos.Y()  << " " << pos.Z() << std::endl;
    std::cout << "direction: " << dir.X() << " " << dir.Y()  << " " << dir.Z() << std::endl;
    std::cout << "momentum:  " << mom << " " << momErr << std::endl;
  }


}

//*************************************************************
void COMET::IGlobalReconModule::FillKinematicsAtTrueVertex(COMET::IHandle<COMET::IG4Trajectory> G4track, 
              COMET::IHandle<COMET::IReconBase> object,
              TLorentzVector& pos, TLorentzVector& posVar, 
              TVector3& dir, TVector3& dirVar, 
              double& mom, double& momErr){
//*************************************************************


  // Propagate the first state in the track to a surface normal to z at the vertex position. 

  COMET::IHandle<COMET::IReconState> tstate = object->GetState(); 
  
  COMET::IHandle<COMET::IReconState> newstate = COMET::IHandle<COMET::IReconState>();

  // The G4 track should exist and the track should have TPC, FGD or P0D components
  if (tstate && G4track && (object->UsesDetector(COMET::IReconBase::kTPC)  || object->UsesDetector(COMET::IReconBase::kFGD) || object->UsesDetector(COMET::IReconBase::kP0D))){    
    bool ok = true;
    TVector3 snorm;
    double dir_max=-100;
    int idir_max=2;
    TVector3 dirS = TrackingUtils::GetDirection(object);
    for (int i=0;i<3;i++){
      if(fabs(dirS(i))> dir_max ){
  dir_max = fabs(dirS(i));
  idir_max=i;
      } 	  
    }
    if (idir_max==0)
      snorm=TVector3(1,0,0);
    else if (idir_max==1)
      snorm=TVector3(0,1,0);
    else if (idir_max==2)
      snorm=TVector3(0,0,1);
    
    if (debug_vertex_extrap)
      std::cout << "Extrapolation to true vertex: dir = (" << dirS.X() << ", " << dirS.Y() << ", " << dirS.Z() << ")" << std::endl;

    if (debug_vertex_extrap)
      COMET::tman("globalAnalysis").SetVerbosity(0,6,0,0);
    
    newstate = COMET::IHandle<COMET::IReconState>(new COMET::ITrackState);
    ok = COMET::tman("globalAnalysis").PropagateToSurface(*tstate,G4track->GetInitialPosition().Vect(),snorm,*newstate);
    
    if (debug_vertex_extrap)
      COMET::tman("globalAnalysis").SetVerbosity(0,0,0,0);
    
  }
  
  FillKinematics(newstate, pos, posVar, dir, dirVar, mom, momErr);

}

//*************************************************************
void COMET::IGlobalReconModule::FillGlobalPIDAlternates(COMET::IHandle<COMET::IG4Trajectory> G4track, COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject){
//*************************************************************


  // __________________________________________________________
  // Fill the Alternates information for inner (P0D+Tracker) Global PIDs.

  if (!(object->UsesDetector(COMET::IReconBase::kTPC)) || object->CheckStatus(COMET::IReconBase::kLikelihoodFit)) return;
  
  // Get the inner object (P0D+TRACKER)
  COMET::IHandle<COMET::IReconBase> innerObject = object;
  if (innerObject->UsesDetector(COMET::IReconBase::kSMRD))
    innerObject = ReconObjectUtils::GetConstituentNotUsingDetector(innerObject,COMET::IReconBase::kSMRD);
  if (innerObject->UsesDetector(COMET::IReconBase::kECal))
    innerObject = ReconObjectUtils::GetConstituentNotUsingDetector(innerObject,COMET::IReconBase::kECal);
  if (innerObject->UsesDetector(COMET::IReconBase::kDSECal))
    innerObject = ReconObjectUtils::GetConstituentNotUsingDetector(innerObject,COMET::IReconBase::kDSECal);
  
  // make sure it is a PID
  COMET::IHandle<COMET::IReconPID> innerPID = innerObject;
  if (!innerPID) return;
  
  

  const COMET::IReconObjectContainer& alter = innerPID->GetAlternates();
  COMET::IReconObjectContainer::const_iterator it;
  for (it=alter.begin();it!=alter.end();it++){  
    // create a new IGlobalPIDAlternate
    IGlobalPIDAlternate* PIDAlternate = new((*globalObject.Alternates)[globalObject.NAlternates++]) IGlobalPIDAlternate;
    FillGlobalPIDAlternate(G4track, *it, *PIDAlternate);
  }

}

//*************************************************************
void COMET::IGlobalReconModule::FillGlobalPIDAlternate(COMET::IHandle<COMET::IG4Trajectory> G4track,
                COMET::IHandle<COMET::IReconBase> object, 
                IGlobalPIDAlternate& PIDAlternate){
//*************************************************************

  // __________________________________________________________
  // Fill the Alternates information for inner (P0D+Tracker) Global PIDs.

  // try to get the state of the object. If it dosen't exists skip the object. 
  try{ COMET::IHandle<COMET::IReconState> tstate = object->GetState();}
  catch(COMET::EReconObject e) { std::cout << "IReconBase with no state. Skip object !!" << std::endl; return; }


  //--------- Information about the detectors that used in this object
  PIDAlternate.Detectors = GetDetectorNumber(object);
  FillDetectorUsed(object, PIDAlternate.DetectorUsed);
    
  //-------  General information 
  
  PIDAlternate.Status     = object->CheckStatus(object->kSuccess);
  PIDAlternate.Chi2       = object->GetQuality();
  PIDAlternate.NDOF       = object->GetNDOF();
  PIDAlternate.Charge     = TrackingUtils::GetCharge(object);	
  PIDAlternate.Length     = ComputeTrackLength(object);

  //------- Front and Back kinematics
  COMET::IHandle<COMET::IReconState> firstState = TrackingUtils::GetFirstState(*object);
  FillKinematics(firstState,
     PIDAlternate.FrontPosition, PIDAlternate.FrontPositionVar,
     PIDAlternate.FrontDirection,PIDAlternate.FrontDirectionVar,
     PIDAlternate.FrontMomentum, PIDAlternate.FrontMomentumError);
  
  COMET::IHandle<COMET::IReconState> lastState  = TrackingUtils::GetLastState(*object);
  FillKinematics(lastState,
     PIDAlternate.BackPosition, PIDAlternate.BackPositionVar,
     PIDAlternate.BackDirection,PIDAlternate.BackDirectionVar,
     PIDAlternate.BackMomentum, PIDAlternate.BackMomentumError);

  //-------- Fill the kinematic extrapolated to the true vertex
  FillKinematicsAtTrueVertex(G4track, object,
           PIDAlternate.PositionAtTrueVertex,  PIDAlternate.PositionVarAtTrueVertex,
           PIDAlternate.DirectionAtTrueVertex, PIDAlternate.DirectionVarAtTrueVertex,
           PIDAlternate.MomentumAtTrueVertex,  PIDAlternate.MomentumErrorAtTrueVertex);
    
  //---------- Set the PID information
  COMET::IHandle<COMET::IReconPID> PID = object;
  if (PID){      
    PIDAlternate.PIDWeight  = PID->GetPIDWeight();
    PIDAlternate.ParticleId = ComputeParticleId(PID);
  }

}

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

  // Default for showers is 0
  double length = 0;
  if (IsTrackLike(object)){
    // Default for tracks
    length = DEFAULT_MAX;
    double length_temp;
    if (object->CheckStatus(COMET::IReconBase::kLikelihoodFit)){ 
      // get the first state and convert it
      COMET::IHandle<COMET::IReconState> firstTState = TrackingUtils::GetFirstState(*object);
      if (!firstTState) return length;
      State state;
      bool ok = COMET::converter().TReconState_to_State(*firstTState,state); 
      if (!ok) return length;

      // get the last node and convert it
      COMET::IHandle<COMET::IReconNode> lastTNode = TrackingUtils::GetLastFittedNode(*object);
      if (!lastTNode) return length;
      Node lastNode;
      ok = COMET::converter().TReconNode_to_Node(lastTNode,lastNode); 
      if (!ok) return length;

      // Propagate from the first to the last node
      ok = COMET::rpman("globalAnalysis").matching_svc().propagate(lastNode,state,length_temp);
      if (ok)
  length = length_temp;
    }
    else{
      Trajectory traj;
      bool ok = COMET::converter().TReconBase_to_Trajectory(object,traj);
      if (ok){ 
  ok = COMET::rpman().matching_svc().compute_length(traj, length_temp);  
  if (ok)
    length = length_temp;
      }
    }
  }

  return length;
}
  
//*************************************************************
void COMET::IGlobalReconModule::FillVertexInfo(COMET::ICOMETEvent& event, COMET::IHandle<COMET::IReconObjectContainer> globalObjects){
//*************************************************************
  //Container names
  TString tname = "ReconTracker";
  TString tvname = "tracker"; // containers are "trakerPrimaryVertex" and "trackerVertices"
  TString gname = "globalVertex";
  TString gvname = "globalRecon";

  //We assume it's global vertices first
  bool global = true;
  TString reconame = gname;
  TString vname = gvname;

  //Fill primary and secondary vertex info
  //**********************************************************
  for(int ibin = 0; ibin < NTIMEBINS; ibin ++){

    TString bin_result_name = reconame + "_timebin";
    bin_result_name += ibin;

    COMET::IHandle<COMET::IAlgorithmResult> ResultBin = event.GetFit(bin_result_name);        
    //Check presence of global vertexing, if not try tracker vertices
    if(ibin == 0 && !ResultBin){
      global = false;
      reconame = tname;
      vname = tvname;
      ResultBin = event.GetFit("ReconTracker_timebin0");
    }

    //If we didn't find a bin, stop looking
    if(!ResultBin) break;

    //Fill tracker to global tracks map
    if(!global){
      TString container_name = "finalTrackerObjects";
      COMET::IHandle<COMET::IReconObjectContainer> finalTracks = ResultBin->GetResultsContainer(container_name);
      if(!finalTracks) continue;
      if(finalTracks->size() == 0) continue;
      fTrackerGlobalIndexMap.clear(); 
      DoAssociationBetweenTrackerAndGlobalObjects(*globalObjects,*finalTracks);
    }

    TString primary_name = vname + "PrimaryVertex";
    COMET::IHandle<COMET::IReconObjectContainer> pvertices = ResultBin->GetResultsContainer(primary_name);
    if(!pvertices) continue;
    if(pvertices->size() == 0) continue;

    COMET::IHandle<COMET::IReconVertex> pv = pvertices->front();
    if(!pv) continue;

    //Fill the primary
    FillVertex(pv, true);

    COMET::IHandle<COMET::IG4PrimaryVertexContainer> g4PrimVert = event.Get<COMET::IG4PrimaryVertexContainer>("truth/G4PrimVertex00");
    if(g4PrimVert) MatchTrueVertex(g4PrimVert);

    TString vertices_name = vname + "Vertices";
    COMET::IHandle<COMET::IReconObjectContainer> svertices =  ResultBin->GetResultsContainer(vertices_name);
    if(!svertices) continue;

    for(COMET::IReconObjectContainer::iterator siter = svertices->begin(); siter != svertices->end(); siter ++){
      
      COMET::IHandle<COMET::IReconVertex> sv = *siter;
      if(!sv) continue;

      //Skip the primary vertex
      if(sv == pv) continue;

      bool vertexCreated = false;

      //Fill secondary vertex
      vertexCreated = FillVertex(sv, false);

      if(g4PrimVert && vertexCreated) MatchTrueVertex(g4PrimVert);
    }
  }  

  fTrackerGlobalIndexMap.clear(); 
}


//*************************************************************
void COMET::IGlobalReconModule::FillGlobalHit(COMET::IHandle<COMET::IHit> hit, IGlobalHit& gHit){
//*************************************************************
  
  if (hit){
    gHit.Position = hit->GetPosition();
    gHit.Time     = hit->GetTime();
    gHit.PositionError = hit->GetUncertainty() ;
    gHit.Charge   = hit->GetCharge();
    gHit.Type     = hit->IsXHit() + (hit->IsYHit())*10 + (hit->IsZHit())*100;
  }  
  else{
    gHit.Position = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    gHit.Time     = DEFAULT_MAX;
    gHit.PositionError = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    gHit.Charge   = DEFAULT_MAX;
    gHit.Type     = DEFAULT_MAX;  
  }

}

//*************************************************************
void COMET::IGlobalReconModule::FillGlobalHit(COMET::IHandle<COMET::IHit> hit, ISMRDHit& smrdHit){
//*************************************************************
  
  if (hit){
    smrdHit.Position = hit->GetPosition();
    smrdHit.Time     = hit->GetTime();
    smrdHit.PositionError = hit->GetUncertainty() ;
    smrdHit.Charge   = hit->GetCharge();
    smrdHit.Type     = hit->IsXHit() + (hit->IsYHit())*10 + (hit->IsZHit())*100;
  }  
  else{
    smrdHit.Position = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    smrdHit.Time     = DEFAULT_MAX;
    smrdHit.PositionError = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    smrdHit.Charge   = DEFAULT_MAX;
    smrdHit.Type     = DEFAULT_MAX; 
  }

}

//*************************************************************
void COMET::IGlobalReconModule::FillSmrdHit(COMET::IHandle<COMET::IHit> hit, ISMRDHit& smrdHit){
//*************************************************************
  enum ESMRDWall {kTop = 0, kBottom, kRight, kLeft};
  //fill general info
  FillGlobalHit(hit, smrdHit);
  if (hit){
    bool isTopWall    = COMET::IGeomInfo::SMRD().IsSMRDTopWall(hit);
    bool isBottomWall = COMET::IGeomInfo::SMRD().IsSMRDBottomWall(hit);
    bool isRightWall  = COMET::IGeomInfo::SMRD().IsSMRDRightWall(hit);
    bool isLeftWall   = COMET::IGeomInfo::SMRD().IsSMRDLeftWall(hit);
    
    if(isTopWall)            
      smrdHit.Wall = kTop;    
    if(isBottomWall)            
      smrdHit.Wall = kBottom;    
    if(isRightWall)            
      smrdHit.Wall = kRight;    
    if(isLeftWall)            
      smrdHit.Wall = kLeft;    		
    smrdHit.Yoke    = COMET::IGeomInfo::SMRD().GetYokeRingNumber(hit);
    smrdHit.Layer   = COMET::IGeomInfo::SMRD().GetSMRDRingNumber(hit);
    smrdHit.Tower   = COMET::IGeomInfo::SMRD().GetSMRDTowerNumber(hit);        
    smrdHit.Counter = COMET::IGeomInfo::SMRD().GetSMRDScintNumber(hit);            	 	
    
  }  
  else{
    //fill with unmeaningful values
    smrdHit.Wall=-1;
    smrdHit.Yoke=-1;
    smrdHit.Layer=-1;
    smrdHit.Tower=-1;
    smrdHit.Counter=-1; 

  }
  
  

}



//*************************************************************
void COMET::IGlobalReconModule::FillSubBaseObject(COMET::IHandle<COMET::IReconBase> object, ISubBaseObject& sub, int det ){
//*************************************************************

  if(ReconObjectUtils::GetOriginalObject(object))
    sub.UniqueID  = ReconObjectUtils::GetOriginalObject(object)->GetUniqueID();
  else
    sub.UniqueID  = object->GetUniqueID();

  sub.Detector  = det;
  sub.Status    = object->CheckStatus(object->kSuccess);
  sub.NDOF      = object->GetNDOF();
  try{ sub.NNodes     = object->GetNodes().size();}
  catch(COMET::EReconObject e) { sub.NNodes = 0; }

  sub.Chi2      = object->GetQuality();
  sub.EDeposit  = TrackingUtils::GetEDeposit(object);


  // The number of constituents. For composite objects
  sub.NConstituents = 0;
  if (object->GetConstituents())
    sub.NConstituents  = object->GetConstituents()->size();

  sub.NHits    = DEFAULT_MAX;
  if (object->GetHits())
    sub.NHits    = object->GetHits()->size();

  // the length of the object for track-like objects
  sub.Length        = ComputeTrackLength(object);

  //--------- Front and back Kinematics
  COMET::IHandle<COMET::IReconState> firstState = TrackingUtils::GetFirstState(*object);
  FillKinematics(firstState,
     sub.FrontPosition, sub.FrontPositionVar,
     sub.FrontDirection,sub.FrontDirectionVar,
     sub.FrontMomentum, sub.FrontMomentumError);
  
  COMET::IHandle<COMET::IReconState> lastState  = TrackingUtils::GetLastState(*object);
  FillKinematics(lastState,
     sub.BackPosition, sub.BackPositionVar,
     sub.BackDirection,sub.BackDirectionVar,
     sub.BackMomentum, sub.BackMomentumError);


  // Fill the true particle when it is found
  double pur, eff;
  COMET::IHandle<COMET::IG4Trajectory> G4track = GetG4Trajectory(*object,pur,eff);
  FillTrueParticle(G4track, pur, eff, sub.TrueParticle);


}

//*************************************************************
void COMET::IGlobalReconModule::FillECALInfo(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject  ){
//*************************************************************


  COMET::IReconBase::Status dets[9] ={COMET::IReconBase::kLeftPECal, COMET::IReconBase::kRightPECal, COMET::IReconBase::kTopPECal, COMET::IReconBase::kBottomPECal, 
           COMET::IReconBase::kLeftTECal, COMET::IReconBase::kRightTECal, COMET::IReconBase::kTopTECal, COMET::IReconBase::kBottomTECal, 
           COMET::IReconBase::kDSECal };

  // Get constituents in any of the calorimeters
  for (int i = 0;i<9;i++){
    if (!object->UsesDetector(dets[i])) continue;
    
    //one needs to take into account that ReconECal doesn`t combine tracks that cross 
    //Left/Right sub-volumes for TECal Top and Bottom 
    COMET::IHandle<COMET::IReconObjectContainer> ecalObjects(new COMET::IReconObjectContainer()); 
   

    ReconObjectUtils::GetConstituentsInDetector(object, dets[i], *ecalObjects);

    for (COMET::IReconObjectContainer::iterator  iter = ecalObjects->begin(); iter!=ecalObjects->end(); iter++){
      COMET::IHandle<COMET::IReconBase> ecalObject = *iter; 
  
      if (!ecalObject) continue;

      IECALObject* ecal; 
      // P0DECAL constituents
      if (i<4){
        if(globalObject.NP0DECALs >= NMAXP0DECAL) continue;          
        ecal = new((*(globalObject.P0DECAL))[globalObject.NP0DECALs++]) IECALObject;
      }
      // TrECAL constituents
      else if (i<8){
        if(globalObject.NECALs >= NMAXECAL) continue;          
        ecal = new((*(globalObject.ECAL))[globalObject.NECALs++]) IECALObject;
        globalObject.NTrECALs++;
      }
      // DsECAL constituents
      else{
        if(globalObject.NECALs >= NMAXECAL) continue;          
        ecal = new((*(globalObject.ECAL))[globalObject.NECALs++]) IECALObject;
        globalObject.NDsECALs++;
      }

      // fill general information
      FillSubBaseObject(ecalObject, *ecal, i+1);
    
      // ECAL specific info
      //      ecal->Cone     = TrackingUtils::GetCone(ecalObject);
      ecal->TrShVal = DEFAULT_MAX;
      ecal->EMHadVal = DEFAULT_MAX;
      ecal->EMEnergy = DEFAULT_MAX;
      ecal->EMEnergyError = DEFAULT_MAX;
      ecal->Cone = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
      ecal->IsShowerLike = false;
      ecal->IsTrackLike = false;      
      ecal->EDeposit = DEFAULT_MAX;
      //      ecal->NHits = 0;
      //      ecal->FrontPosition = TLorentzVector(0,0,0,0);
      //      ecal->FrontDirection = TVector3(0,0,0);
    
      // Perkin 20110823
      // ECal layer info (from ITrackerReconECalModule.cxx)
      // Defalut values -1
      ecal->mostUpStreamLayerHit   = -1;
      ecal->mostDownStreamLayerHit = -1;
    
    
      COMET::IHandle<COMET::IReconObjectContainer> conContainer = ecalObject->GetConstituents();
      if (conContainer){
        // The  PID variables are not saved in the main object but in its constituent
        COMET::IHandle<COMET::IReconBase> ecalObject2 = *(conContainer->begin());
      
        if (ecalObject2){ 
  
    std::string constituent_class_name =  ecalObject2->ClassName();
  
    // Looking for Pids with a Constituent of a Track first.
    if(constituent_class_name.find("COMET::IReconTrack") != std::string::npos) {
      COMET::IHandle<COMET::IReconTrack> track = ecalObject2;
    
      if(track->GetEDeposit()) {
        // These bools are to show what the ReconECal using default 
        // cuts thinks the object is, a shower, or a track.
        ecal->IsShowerLike = false;
        ecal->IsTrackLike  = true;
      
        ecal->EDeposit = track->GetEDeposit();
        ecal->Width    = track->GetWidth();
        //	      ecal->FrontPosition = track->GetPosition();
        //	      ecal->FrontDirection = track->GetDirection();
      
        // Perkin 20110823
        // ECal layer info (from ITrackerReconECalModule.cxx)
      
        if(track->GetHits()){
          ecal->mostUpStreamLayerHit   = 100;
          ecal->mostDownStreamLayerHit = -1;
        
          COMET::IHandle<COMET::IHitSelection> hitsel = track->GetHits();
          for(COMET::IHitSelection::iterator hit = hitsel->begin(); 
        hit != hitsel->end(); ++hit){
    
      COMET::IHandle<COMET::IReconHit> reconHit = (*hit);
      TVector3 HitPos = reconHit->GetContributor(0)->GetPosition();
      Int_t layerHit = COMET::IGeomInfo::Get().ECAL().GetLayerNumber(HitPos);
    
      if(layerHit < ecal->mostUpStreamLayerHit)
        ecal->mostUpStreamLayerHit = layerHit;
    
      if(layerHit > ecal->mostDownStreamLayerHit)
        ecal->mostDownStreamLayerHit = layerHit;
          }
        }
      }
    
    }
    else if (constituent_class_name.find("COMET::IReconShower") != std::string::npos) {
      COMET::IHandle<COMET::IReconShower> shower = ecalObject2;
    
      if(shower->GetEDeposit()) {
        ecal->IsShowerLike = true;
        ecal->IsTrackLike = false;
      
        ecal->EDeposit = shower->GetEDeposit();
        ecal->Cone = shower->GetConeAngle();

        //	      ecal->FrontPosition = shower->GetPosition();
        //              ecal->FrontDirection = shower->GetDirection();
      
        // Perkin 20110823
        // ECal layer info (from ITrackerReconECalModule.cxx)
        if(shower->GetHits()){
          ecal->mostUpStreamLayerHit   = 100;
          ecal->mostDownStreamLayerHit = -1;
      
          std::vector<TVector3> layerHitPos(36);
          std::vector<double> layerCharge(36,0.0);
        

          COMET::IHandle<COMET::IHitSelection> hitsel = shower->GetHits();
          for(COMET::IHitSelection::iterator hit = hitsel->begin(); 
        hit != hitsel->end(); ++hit){
    
      COMET::IHandle<COMET::IReconHit> reconHit = (*hit);
      TVector3 HitPos = reconHit->GetContributor(0)->GetPosition();
      Int_t layerHit = COMET::IGeomInfo::Get().ECAL().GetLayerNumber(HitPos);
    
      if(layerHit < ecal->mostUpStreamLayerHit)
        ecal->mostUpStreamLayerHit = layerHit;
    
      if(layerHit > ecal->mostDownStreamLayerHit)
        ecal->mostDownStreamLayerHit = layerHit;

      // Charge weighted layer positon
      TVector3 chargePos = reconHit->GetPosition() * reconHit->GetCharge();
      layerHitPos[layerHit] += chargePos;
      layerCharge[layerHit] += reconHit->GetCharge();
          }

          if ( ecal->mostDownStreamLayerHit!=-1 &&  ecal->mostDownStreamLayerHit!=-1){
    
      // Charge weighted layer position
      for(int i = 0; i < 36; ++i){
        if(layerCharge[i] != 0.0){
          layerHitPos[i] *= (1 / layerCharge[i]);
        }
      }
    
      // For showers the length is the distance between the the first and last layer hits
      TVector3 objLength = layerHitPos[ecal->mostDownStreamLayerHit] - layerHitPos[ecal->mostUpStreamLayerHit];
      ecal->Length = objLength.Mag();
                  if(ecalObject2->Get< COMET::IRealDatum >("ShowerAngle")){
                      ecal->Width.SetX(ecal->Length*tan(ecalObject2->Get< COMET::IRealDatum >("ShowerAngle")->GetValue()));
                  }
                  else{
                      ecal->Width.SetX(DEFAULT_MAX);
                  }
          }   
        }
      }    
      
      }
  
  
    if(ecalObject2->Get< COMET::IRealDatum >("TrShval")){ 
      ecal->TrShVal = ecalObject2->Get< COMET::IRealDatum >("TrShval")->GetValue();
    }
    else{
      ecal->TrShVal = -1;
    }
  
    if(ecalObject2->Get< COMET::IRealDatum >("EMHadVal")){ 
      ecal->EMHadVal = ecalObject2->Get< COMET::IRealDatum >("EMHadVal")->GetValue();
    }
    else{
      ecal->EMHadVal = -1;
    }
  
    // Get the EM Energy Fit Result and Uncertainty
    COMET::IHandle<COMET::IReconObjectContainer> conContainer2 = ecalObject2->GetConstituents();
    if (conContainer2){ 
      COMET::IHandle<COMET::IReconBase> ecalObject3 = *(conContainer2->begin());
      if (ecalObject3){ 
        std::string cluster_constituent_class_name =  ecalObject3->ClassName();
        if(cluster_constituent_class_name.find("COMET::IReconCluster") != std::string::npos) {
          COMET::IHandle<COMET::IReconCluster> cluster = ecalObject3;
        
          double totalHitCharge = 0;
          // Get the number of hits for the ecal cluster
          for (COMET::IHitSelection::const_iterator it = cluster->GetHits()->begin(); it != cluster->GetHits()->end(); ++it) {
      ecal->AverageHitTime += (*it)->GetTime() * (*it)->GetCharge();
      totalHitCharge += (*it)->GetCharge();
          }
          ecal->AverageHitTime /= totalHitCharge;
        
          // EM Energy Fit Result   
          if(cluster->Get< COMET::IRealDatum >("EMEnergyFitResult")){ 
      ecal->EMEnergy = cluster->Get< COMET::IRealDatum >("EMEnergyFitResult")->GetVector().front();
          }
          else{
      ecal->EMEnergy = -1;
          }
        
          // EM Energy Fit Uncertainty
          if(cluster->Get< COMET::IRealDatum >("EMEnergyFitUncertainty")){ 
      ecal->EMEnergyError = cluster->Get< COMET::IRealDatum >("EMEnergyFitUncertainty")->GetValue();
          }
          else{
      ecal->EMEnergyError = -1;
          }
          //
          // EM Energy Fit Likelihoods
          if(cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")){
      ecal->EMEnergy_Likelihood_energyGrad = cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")->GetVector()[0];
      ecal->EMEnergy_Likelihood_qsum_like = cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")->GetVector()[1];
      ecal->EMEnergy_Likelihood_energy_qsumGrad = cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")->GetVector()[2];
          }
          else{
      ecal->EMEnergy_Likelihood_energyGrad = -1;
      ecal->EMEnergy_Likelihood_energy_qsumGrad = -1;
      ecal->EMEnergy_Likelihood_qsum_like = -1;
          }
        
          // EM Energy Fit Parameters
          if(cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")){
      ecal->EMEnergy_Para_QSum = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[0];
      ecal->EMEnergy_Para_QMean = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[1];
      ecal->EMEnergy_Para_QRMS = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[2];
      ecal->EMEnergy_Para_QSkew = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[3];
      ecal->EMEnergy_Para_QMax = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[4];
          }
          else{
      ecal->EMEnergy_Para_QSum = -1;
      ecal->EMEnergy_Para_QMean = -1;
      ecal->EMEnergy_Para_QRMS = -1;
      ecal->EMEnergy_Para_QSkew = -1;
      ecal->EMEnergy_Para_QMax = -1;
          }
          //
          // PID output from Kalman Fitter
          if(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")){
      ecal->KFParameter     = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[0];
      //ecal->KFMultiTracksTPC  = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[1];
      //ecal->KFMultiTracksECAL = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[2];
      ecal->KFNNodes = (int)cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[3];
      ecal->KFParameterNodes = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[4];
      ecal->KFNDOF = (int)cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[5];
      ecal->KFQuality = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[6];
      //ecal->KFIsMatched = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[7];
      ecal->KFHitQuality = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[8];
          }
          else{
      ecal->KFParameter     = -1;
      //ecal->KFMultiTracksTPC  = -1;
      //ecal->KFMultiTracksECAL = -1;
      ecal->KFNNodes = -1;
      ecal->KFParameterNodes = -1;
      ecal->KFQuality = -1;
      ecal->KFNDOF = -1;
      //ecal->KFIsMatched = false;
      ecal->KFHitQuality = -1;
          }
          //
          // Michel Electron Tag Output
          if(cluster->Get< COMET::IRealDatum >("NDelayedClusters")){
      ecal->MElectronTag_NDelayedCluster = cluster->Get< COMET::IRealDatum >("NDelayedClusters")->GetVector()[0];
          }
          else{
      ecal->MElectronTag_NDelayedCluster= -1;
          } 
        
          if(ecal->MElectronTag_NDelayedCluster > 0){
      for(Int_t dclust = 0; dclust < ecal->MElectronTag_NDelayedCluster; ++dclust){
      
        std::string delayed_name = "DelayedCluster_";
        char dclust_string [3];
        sprintf(dclust_string, "%d", dclust);
        std::string delayed_clust = (delayed_name + dclust_string);
      
        if(cluster->Get<COMET::IRealDatum>(delayed_clust.c_str())){
        
          std::vector<double> Delta = cluster->Get<COMET::IRealDatum>(delayed_clust.c_str())->GetVector();
        
          //Grab Michel Candidate Variabes from Datum
          ecal->MElectronTag_NDelayedHits.push_back(Delta[0]);
          ecal->MElectronTag_DeltaX.push_back(Delta[1]);
          ecal->MElectronTag_DeltaY.push_back(Delta[2]);
          ecal->MElectronTag_DeltaZ.push_back(Delta[3]);
          ecal->MElectronTag_DeltaT.push_back(Delta[4]);
          ecal->MElectronTag_EDep.push_back(Delta[5]);
          ecal->MElectronTag_TrackEnd.push_back(Delta[6]);
        }
      }
          }// Michel Tag
        }
      }
    
      // Leigh Whitehead 06/09/10 
      // Add in the PID neural network input variables.
    
      // Get the AMR
      if(ecalObject2->Get< COMET::IRealDatum >("AMR")) {
        ecal->PID_AMR = ecalObject2->Get< COMET::IRealDatum >("AMR")->GetValue();
      }
      else{
        ecal->PID_AMR = -1;
      }
    
      // Get the Angle
      if(ecalObject2->Get< COMET::IRealDatum >("Angle")){ 
        ecal->PID_Angle = ecalObject2->Get< COMET::IRealDatum >("Angle")->GetValue();
      }
      else{
        ecal->PID_Angle = -1;
      }
    
      // Get the Max Ratio
      if(ecalObject2->Get< COMET::IRealDatum >("Max_Ratio")){ 
        ecal->PID_Max_Ratio = ecalObject2->Get< COMET::IRealDatum >("Max_Ratio")->GetValue();
      }
      else {
        ecal->PID_Max_Ratio = -1;
      }
    
      // Get the Shower Angle
      if(ecalObject2->Get< COMET::IRealDatum >("ShowerAngle")){ 
        ecal->PID_ShowerAngle = ecalObject2->Get< COMET::IRealDatum >("ShowerAngle")->GetValue();
      }
      else{
        ecal->PID_ShowerAngle = -1;
      }
    
      // Asymmetry
      if(ecalObject2->Get< COMET::IRealDatum >("Asymmetry")){
        ecal->PID_Asymmetry = ecalObject2->Get< COMET::IRealDatum >("Asymmetry")->GetValue();
      }
      else{
        ecal->PID_Asymmetry = -1;
      }
    
      // Mean Position
      if(ecalObject2->Get< COMET::IRealDatum >("MeanPosition")) {
        ecal->PID_MeanPosition = ecalObject2->Get< COMET::IRealDatum >("MeanPosition")->GetValue();
      }
      else{
        ecal->PID_MeanPosition = -1;
      }
    
      // Qskew - same info also saved in the QSkew variable for the 
      // EMEnergyFit. No need to save it here as well
      /*
      if(ecalObject2->Get< COMET::IRealDatum >("Qskew")){ 
        ecal->PID_Qskew = ecalObject2->Get< COMET::IRealDatum >("Qskew")->GetValue();
      }
      else {
        ecal->PID_Qskew = -1;
      }
      */
      // Shower Width
      if(ecalObject2->Get< COMET::IRealDatum >("ShowerWidth")){ 
        ecal->PID_ShowerWidth = ecalObject2->Get< COMET::IRealDatum >("ShowerWidth")->GetValue();
        ecal->Width.SetX(ecal->PID_ShowerWidth);
      }
      else{
        ecal->PID_ShowerWidth = -1;
        ecal->Width = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
      }
    
      // EM Likelihood  - same info also saved for the 
      // EMEnergyFit. No need to save it here as well
      /*
      if(ecalObject2->Get< COMET::IRealDatum >("EMLikelihood")){ 
        ecal->PID_EMLikelihood = ecalObject2->Get< COMET::IRealDatum >("EMLikelihood")->GetValue();
      }
      else{
        ecal->PID_EMLikelihood = -1;
      }
      */
    
    }	 
   
        }
      }
    }
  }
 }

//*************************************************************
void COMET::IGlobalReconModule::FillTrackerInfo(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject){
//*************************************************************

  // We created dummy objects for the 18 possible Tracker combinations:


  // No detector missing
  // TPC1-FGD1                
  // TPC1-FGD1-TPC2
  // TPC1-FGD1-TPC2-FGD2
  // TPC1-FGD1-TPC2-FGD2-TPC3
  // FGD1-TPC2
  // FGD1-TPC2-FGD2
  // FGD1-TPC2-FGD2-TPC3
  // TPC2-FGD2
  // TPC2-FGD2-TPC3
  // FGD2-TPC3
 
  COMET::IReconTrack dummy[10];
  dummy[0].AddDetector(COMET::IReconBase::kTPC1);
  dummy[0].AddDetector(COMET::IReconBase::kFGD1);

  dummy[1].AddDetector(COMET::IReconBase::kTPC1);
  dummy[1].AddDetector(COMET::IReconBase::kFGD1);
  dummy[1].AddDetector(COMET::IReconBase::kTPC2);

  dummy[2].AddDetector(COMET::IReconBase::kTPC1);
  dummy[2].AddDetector(COMET::IReconBase::kFGD1);
  dummy[2].AddDetector(COMET::IReconBase::kTPC2);
  dummy[2].AddDetector(COMET::IReconBase::kFGD2);

  dummy[3].AddDetector(COMET::IReconBase::kTPC1);
  dummy[3].AddDetector(COMET::IReconBase::kFGD1);
  dummy[3].AddDetector(COMET::IReconBase::kTPC2);
  dummy[3].AddDetector(COMET::IReconBase::kFGD2);
  dummy[3].AddDetector(COMET::IReconBase::kTPC3);

  dummy[4].AddDetector(COMET::IReconBase::kFGD1);
  dummy[4].AddDetector(COMET::IReconBase::kTPC2);

  dummy[5].AddDetector(COMET::IReconBase::kFGD1);
  dummy[5].AddDetector(COMET::IReconBase::kTPC2);
  dummy[5].AddDetector(COMET::IReconBase::kFGD2);

  dummy[6].AddDetector(COMET::IReconBase::kFGD1);
  dummy[6].AddDetector(COMET::IReconBase::kTPC2);
  dummy[6].AddDetector(COMET::IReconBase::kFGD2);
  dummy[6].AddDetector(COMET::IReconBase::kTPC3);

  dummy[7].AddDetector(COMET::IReconBase::kTPC2);
  dummy[7].AddDetector(COMET::IReconBase::kFGD2);

  dummy[8].AddDetector(COMET::IReconBase::kTPC2);
  dummy[8].AddDetector(COMET::IReconBase::kFGD2);
  dummy[8].AddDetector(COMET::IReconBase::kTPC3);

  dummy[9].AddDetector(COMET::IReconBase::kFGD2);
  dummy[9].AddDetector(COMET::IReconBase::kTPC3);


  // one FGD missing
  // TPC1-FGD1-TPC2-TPC3
  // TPC1-FGD2-TPC2-TPC3
  // TPC1-TPC2
  // TPC2-TPC3

  dummy[10].AddDetector(COMET::IReconBase::kTPC1);
  dummy[10].AddDetector(COMET::IReconBase::kFGD1);
  dummy[10].AddDetector(COMET::IReconBase::kTPC2);
  dummy[10].AddDetector(COMET::IReconBase::kTPC3);

  dummy[11].AddDetector(COMET::IReconBase::kTPC1);
  dummy[11].AddDetector(COMET::IReconBase::kTPC2);
  dummy[11].AddDetector(COMET::IReconBase::kFGD2);
  dummy[11].AddDetector(COMET::IReconBase::kTPC3);

  dummy[12].AddDetector(COMET::IReconBase::kTPC1);
  dummy[12].AddDetector(COMET::IReconBase::kTPC2);

  dummy[13].AddDetector(COMET::IReconBase::kTPC2);
  dummy[13].AddDetector(COMET::IReconBase::kTPC3);


  // two FGDs missing
  // TPC1-TPC2-TPC3

  dummy[14].AddDetector(COMET::IReconBase::kTPC1);
  dummy[14].AddDetector(COMET::IReconBase::kTPC2);
  dummy[14].AddDetector(COMET::IReconBase::kTPC3);


  // one TPC missing
  // TPC1-FGD1-TPC3
  // TPC1-FGD2-TPC3

  dummy[15].AddDetector(COMET::IReconBase::kTPC1);
  dummy[15].AddDetector(COMET::IReconBase::kFGD1);
  dummy[15].AddDetector(COMET::IReconBase::kTPC3);

  dummy[16].AddDetector(COMET::IReconBase::kTPC1);
  dummy[16].AddDetector(COMET::IReconBase::kFGD2);
  dummy[16].AddDetector(COMET::IReconBase::kTPC3);


  // two FGDs and one TPC missing
  // TPC1-TPC3

  dummy[17].AddDetector(COMET::IReconBase::kTPC1);
  dummy[17].AddDetector(COMET::IReconBase::kTPC3);



  for (int i = 0;i<18;i++){
    if (!object->UsesDetector(dummy[i].GetDetectors())) continue;

    COMET::IHandle<COMET::IReconBase> trackerObject = ReconObjectUtils::GetConstituentInDetector(object, dummy[i].GetDetectors());

    if (trackerObject && globalObject.NTRACKERs < NMAXTRACKER){    
      
      ITrackerObject* tracker = new((*(globalObject.TRACKER))[globalObject.NTRACKERs++]) ITrackerObject;

      // fill general information
      int det = GetDetectorNumber(trackerObject);
      FillSubBaseObject(trackerObject, *tracker,det);

      // Fill charge for object.
      tracker->Charge = TrackingUtils::GetCharge(trackerObject);	

    }
  }
}

//*************************************************************
void COMET::IGlobalReconModule::FillTPCInfo(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject  ){
//*************************************************************

  COMET::IReconBase::Status dets[3] ={COMET::IReconBase::kTPC1, COMET::IReconBase::kTPC2, COMET::IReconBase::kTPC3};
  
  for (int i = 0;i<3;i++){
    if (!object->UsesDetector(dets[i])) continue;
    
    COMET::IHandle<COMET::IReconBase> tpcObject = ReconObjectUtils::GetConstituentInDetector(object, dets[i]);

    if (tpcObject && globalObject.NTPCs < NMAXTPC){    
      
      ITPCObject* tpc = new((*(globalObject.TPC))[globalObject.NTPCs++]) ITPCObject;
      FillTPCObject(tpcObject, *tpc,i+1);

    }
  }
}

//*************************************************************
void COMET::IGlobalReconModule::FillTPCObject(COMET::IHandle<COMET::IReconBase> tpcObject, ITPCObject& tpc, int det){
//*************************************************************

  // fill general information
  FillSubBaseObject(tpcObject, tpc, det);
  
  tpc.Charge        = TrackingUtils::GetCharge(tpcObject);	
  
  tpc.NTrun=0;
  tpc.Ccorr=DEFAULT_MAX;
  tpc.PullEle=DEFAULT_MAX;
  tpc.PullMuon=DEFAULT_MAX;
  tpc.PullPion=DEFAULT_MAX;
  tpc.PullKaon=DEFAULT_MAX;
  tpc.PullProton=DEFAULT_MAX;
  
  tpc.SigmaEle=DEFAULT_MAX;
  tpc.SigmaMuon=DEFAULT_MAX;
  tpc.SigmaPion=DEFAULT_MAX;
  tpc.SigmaKaon=DEFAULT_MAX;
  tpc.SigmaProton=DEFAULT_MAX;
  
  tpc.dEdxexpEle=DEFAULT_MAX;
  tpc.dEdxexpMuon=DEFAULT_MAX;
  tpc.dEdxexpPion=DEFAULT_MAX;
  tpc.dEdxexpKaon=DEFAULT_MAX;
  tpc.dEdxexpProton=DEFAULT_MAX;
  
  // fill PID information when exists 
  COMET::IHandle<COMET::IReconPID> pid = tpcObject;
  
  if (pid){    
    bool isKalman = pid->CheckStatus(COMET::IReconBase::kKalmanFit);
    if (isKalman){
       COMET::IHandle<COMET::IReconObjectContainer> LikFits = pid->GetConstituents();
       for( COMET::IReconObjectContainer::iterator  lk = LikFits->begin(); lk != LikFits->end();lk++) {
          COMET::IHandle<COMET::IReconPID> lkPiD = *lk;
          COMET::IHandle<COMET::IReconTrack> lkTrack = *(lkPiD->Get<COMET::IReconObjectContainer>("constituents")->begin());
          if( (lkTrack)->Get< COMET::IRealDatum >("NHorRows") && (lkTrack)->Get< COMET::IRealDatum >("NVerRows") ){
      tpc.NVerRows.push_back(int((lkTrack)->Get< COMET::IRealDatum >("NVerRows")->GetValue()));
      tpc.NHorRows.push_back(int((lkTrack)->Get< COMET::IRealDatum >("NHorRows")->GetValue()));
          }
          if ( (lkTrack)->Get< COMET::IRealDatum >("T0Source") ){
      tpc.T0Source.push_back(int((lkTrack)->Get< COMET::IRealDatum >("T0Source")->GetValue()));
          }
       }
    } else{    
       COMET::IHandle<COMET::IReconTrack> trackTrack = *(pid->Get<COMET::IReconObjectContainer>("constituents")->begin());
       if( (trackTrack)->Get< COMET::IRealDatum >("NHorRows") && (trackTrack)->Get< COMET::IRealDatum >("NVerRows") ){
   tpc.NVerRows.push_back(int((trackTrack)->Get< COMET::IRealDatum >("NVerRows")->GetValue()));
   tpc.NHorRows.push_back(int((trackTrack)->Get< COMET::IRealDatum >("NHorRows")->GetValue()));
       }
       if ( (trackTrack)->Get< COMET::IRealDatum >("T0Source") ){
         tpc.T0Source.push_back(int((trackTrack)->Get< COMET::IRealDatum >("T0Source")->GetValue()));
       }
    }

    // Changed from tpcPid_PullEle to NTrun to fix bug 569 - no reason why this should work, 
    // but it did, so the bug still needs investigating.

    COMET::IHandle<COMET::IRealDatum> datum_ele=(pid)->Get<COMET::IRealDatum>("tpcPid_NTrun");
    if(datum_ele){
      tpc.NTrun         = pid->Get<COMET::IRealDatum>("tpcPid_NTrun")->GetValue();
      tpc.Ccorr         = pid->Get<COMET::IRealDatum>("tpcPid_Ccorr")->GetValue();
      tpc.PullEle       = pid->Get<COMET::IRealDatum>("tpcPid_PullEle")->GetValue();
      tpc.PullMuon      = pid->Get<COMET::IRealDatum>("tpcPid_PullMuon")->GetValue();
      tpc.PullPion      = pid->Get<COMET::IRealDatum>("tpcPid_PullPion")->GetValue();
      tpc.PullKaon      = pid->Get<COMET::IRealDatum>("tpcPid_PullKaon")->GetValue();
      tpc.PullProton    = pid->Get<COMET::IRealDatum>("tpcPid_PullProton")->GetValue();
      tpc.dEdxexpEle    = pid->Get<COMET::IRealDatum>("tpcPid_dEdxexpEle")->GetValue();
      tpc.dEdxexpMuon   = pid->Get<COMET::IRealDatum>("tpcPid_dEdxexpMuon")->GetValue();
      tpc.dEdxexpPion   = pid->Get<COMET::IRealDatum>("tpcPid_dEdxexpPion")->GetValue();
      tpc.dEdxexpKaon   = pid->Get<COMET::IRealDatum>("tpcPid_dEdxexpKaon")->GetValue();
      tpc.dEdxexpProton = pid->Get<COMET::IRealDatum>("tpcPid_dEdxexpProton")->GetValue();
      tpc.SigmaEle      = pid->Get<COMET::IRealDatum>("tpcPid_SigmaEle")->GetValue();
      tpc.SigmaMuon     = pid->Get<COMET::IRealDatum>("tpcPid_SigmaMuon")->GetValue();
      tpc.SigmaPion     = pid->Get<COMET::IRealDatum>("tpcPid_SigmaPion")->GetValue();
      tpc.SigmaKaon     = pid->Get<COMET::IRealDatum>("tpcPid_SigmaKaon")->GetValue();
      tpc.SigmaProton   = pid->Get<COMET::IRealDatum>("tpcPid_SigmaProton")->GetValue();
    }
  }

}


//*************************************************************
void COMET::IGlobalReconModule::FillFGDInfo(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject  ){
//*************************************************************

  // Get the tracker recon version of the this object, if it is a FGD iso-recon object.
  // The global recon version of the object lacks the PID information.
  COMET::IHandle<COMET::IReconBase> newFgdIso = GetTrackerReconVersionOfFGDIsoTrack(object); 	 
  if (newFgdIso) { 	 
    object = newFgdIso; 	 
  }

  COMET::IReconBase::Status dets[2] ={COMET::IReconBase::kFGD1, COMET::IReconBase::kFGD2};

  for (int i = 0;i<2;i++){
    if (!object->UsesDetector(dets[i])) continue;

    COMET::IHandle<COMET::IReconBase> fgdObject =ReconObjectUtils::GetConstituentInDetector(object, dets[i]);

    if (fgdObject && globalObject.NFGDs < NMAXFGD){    

      IFGDObject* fgd = new((*(globalObject.FGD))[globalObject.NFGDs++]) IFGDObject;

      // fill general information
      FillSubBaseObject(fgdObject, *fgd,i+1);

      // fill average hit time
      COMET::IHandle<COMET::IRealDatum> fgdTimeRD = fgdObject->Get<COMET::IRealDatum>("averageHitTime");
      if (fgdTimeRD) {
  fgd->avgtime = fgdTimeRD->GetValue();
      } else {
  fgd->avgtime = DEFAULT_MAX;
      }

      // Fill charge per layer information
      for (int ilayer=0; ilayer<30; ilayer++) {
  fgd->chargePerLayer[ilayer] = 0.;
  fgd->chargePerLayerAttenCorr[ilayer] = 0.;
      }
      std::map<int,double> chargeMap;
      if (fgdUtils::GetChargeByLayers(fgdObject->GetHits(),chargeMap)) {
  for (std::map<int,double>::iterator it=chargeMap.begin() ; it != chargeMap.end(); it++ ) {
    int layer = (*it).first;
    if (i==1) {
      layer -= 30;
    }
    fgd->chargePerLayer[layer] = (*it).second;
  }
      }
      std::map<int,double> chargeMapAC;
      if (fgdUtils::GetFibreCorrectedChargeByLayersReconFit(fgdObject->GetHits(),fgdObject,chargeMapAC)) {
  for (std::map<int,double>::iterator it=chargeMapAC.begin() ; it != chargeMapAC.end(); it++ ) {
    int layer = (*it).first;
    if (i==1) {
      layer -= 30;
    }
    fgd->chargePerLayerAttenCorr[layer] = (*it).second;
  }
      }

//      for (int ilayer=0; ilayer<30; ilayer++) {
//	std::cout << "track in fgd " << i+1 << "; charge per layer[" << ilayer << "]: (raw,attencorr) = (" << fgd->chargePerLayer[ilayer] << "," << fgd->chargePerLayerAttenCorr[ilayer] << ")" << std::endl;
//      }

      // fill PID Information

      // fully contained info
      COMET::IHandle<COMET::IIntegerDatum> RD_geomFC = fgdObject->Get< COMET::IIntegerDatum >("fgdPidPar_geomFC");
      fgd->fgdContainment  = 0;
      if (RD_geomFC) {
  fgd->fgdContainment  = RD_geomFC->GetValue();
      }

      // measured E and x
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_E = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_E");
      fgd->E = -1.;
      if (RD_fgdPidPar_E) {
  fgd->E = RD_fgdPidPar_E->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_x = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_x");
      fgd->x = -1.;
      if (RD_fgdPidPar_x) {
  fgd->x = RD_fgdPidPar_x->GetValue();
      }

      // expected E at the measured x
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_E_exp_EvsXPull_Muon = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_E_exp_EvsXPull_Muon");
      fgd->E_exp_muon = -1.;
      if (RD_fgdPidPar_E_exp_EvsXPull_Muon) {
  fgd->E_exp_muon = RD_fgdPidPar_E_exp_EvsXPull_Muon->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_E_exp_EvsXPull_Pion = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_E_exp_EvsXPull_Pion");
      fgd->E_exp_pion = -1.;
      if (RD_fgdPidPar_E_exp_EvsXPull_Pion) {
  fgd->E_exp_pion = RD_fgdPidPar_E_exp_EvsXPull_Pion->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_E_exp_EvsXPull_Proton = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_E_exp_EvsXPull_Proton");
      fgd->E_exp_proton = -1.;
      if (RD_fgdPidPar_E_exp_EvsXPull_Proton) {
  fgd->E_exp_proton = RD_fgdPidPar_E_exp_EvsXPull_Proton->GetValue();
      }

      // sigma E at the measured x
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_error_EvsXPull_Muon = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_error_EvsXPull_Muon");
      fgd->sigmaE_muon = -1.;
      if (RD_fgdPidPar_error_EvsXPull_Muon) {
  fgd->sigmaE_muon = RD_fgdPidPar_error_EvsXPull_Muon->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_error_EvsXPull_Pion = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_error_EvsXPull_Pion");
      fgd->sigmaE_pion = -1.;
      if (RD_fgdPidPar_error_EvsXPull_Pion) {
  fgd->sigmaE_pion = RD_fgdPidPar_error_EvsXPull_Pion->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidPar_error_EvsXPull_Proton = fgdObject->Get< COMET::IRealDatum >("fgdPidPar_error_EvsXPull_Proton");
      fgd->sigmaE_proton = -1.;
      if (RD_fgdPidPar_error_EvsXPull_Proton) {
  fgd->sigmaE_proton = RD_fgdPidPar_error_EvsXPull_Proton->GetValue();
      }

      // Pull values
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidWgt_EvsXPull_Muon = fgdObject->Get< COMET::IRealDatum >("fgdPidWgt_EvsXPull_Muon");
      fgd->PullMuon = -1.e10;
      if (RD_fgdPidWgt_EvsXPull_Muon) {
  fgd->PullMuon = RD_fgdPidWgt_EvsXPull_Muon->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidWgt_EvsXPull_Pion = fgdObject->Get< COMET::IRealDatum >("fgdPidWgt_EvsXPull_Pion");
      fgd->PullPion = -1.e10;
      if (RD_fgdPidWgt_EvsXPull_Pion) {
  fgd->PullPion = RD_fgdPidWgt_EvsXPull_Pion->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidWgt_EvsXPull_Proton = fgdObject->Get< COMET::IRealDatum >("fgdPidWgt_EvsXPull_Proton");
      fgd->PullProton = -1.e10;
      if (RD_fgdPidWgt_EvsXPull_Proton) {
  fgd->PullProton = RD_fgdPidWgt_EvsXPull_Proton->GetValue();
      }
      COMET::IHandle<COMET::IRealDatum> RD_fgdPidWgt_EvsXPull_NotSet = fgdObject->Get< COMET::IRealDatum >("fgdPidWgt_EvsXPull_NotSet");
      fgd->PullNotSet = -1.e10;
      if (RD_fgdPidWgt_EvsXPull_NotSet) {
  fgd->PullNotSet = RD_fgdPidWgt_EvsXPull_NotSet->GetValue();
      }

     //Vertex activity values
      int fgdNum = 2;
      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_FgdN;
      RD_fgdVA_FgdN = object->Get<COMET::IRealDatum>("fgdVA_FgdN");
      if( RD_fgdVA_FgdN ){
   fgdNum = int( RD_fgdVA_FgdN->GetValue() );
      }

      fgd->isFgdVA = 0;
      if( fgdNum == i )
        fgd->isFgdVA = 1;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_flag;
      RD_fgdVA_flag = object->Get<COMET::IRealDatum>("fgdVA_flag");
      if( RD_fgdVA_flag && i == fgdNum ){
   fgd->isFgdVA_flag = int( RD_fgdVA_flag->GetValue() );
      }
      else
        fgd->isFgdVA_flag = 0;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_upMinZ;
      RD_fgdVA_upMinZ = object->Get<COMET::IRealDatum>("fgdVA_upMinZ");
      if( RD_fgdVA_upMinZ && i == fgdNum ){
   fgd->fgdVA_upMinZ = RD_fgdVA_upMinZ->GetValue();
      }
      else
        fgd->fgdVA_upMinZ = 0;
      
      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_downMaxZ;
      RD_fgdVA_downMaxZ = object->Get<COMET::IRealDatum>("fgdVA_downMaxZ");
      if( RD_fgdVA_downMaxZ && i == fgdNum ){
   fgd->fgdVA_downMaxZ = RD_fgdVA_downMaxZ->GetValue();
      }
      else
        fgd->fgdVA_downMaxZ = 0;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_upCosTheta;
      RD_fgdVA_upCosTheta = object->Get<COMET::IRealDatum>("fgdVA_upCosTheta");
      if( RD_fgdVA_upCosTheta && i == fgdNum ){
   fgd->fgdVA_upCosTheta = RD_fgdVA_upCosTheta->GetValue();
      }
      else
        fgd->fgdVA_upCosTheta = 0;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_downCosTheta;
      RD_fgdVA_downCosTheta = object->Get<COMET::IRealDatum>("fgdVA_downCosTheta");
      if( RD_fgdVA_downCosTheta && i == fgdNum ){
   fgd->fgdVA_downCosTheta = RD_fgdVA_downCosTheta->GetValue();
      }
      else
        fgd->fgdVA_downCosTheta = 0;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_otherUpQ;
      RD_fgdVA_otherUpQ = object->Get<COMET::IRealDatum>("fgdVA_otherUpQ");
      if( RD_fgdVA_otherUpQ && i == fgdNum ){
   fgd->fgdVA_otherUpQ = RD_fgdVA_otherUpQ->GetValue();
      }
      else
        fgd->fgdVA_otherUpQ = 0;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_otherDownQ;
      RD_fgdVA_otherDownQ = object->Get<COMET::IRealDatum>("fgdVA_otherDownQ");
      if( RD_fgdVA_otherDownQ && i == fgdNum ){
   fgd->fgdVA_otherDownQ = RD_fgdVA_otherDownQ->GetValue();
      }
      else
        fgd->fgdVA_otherDownQ = 0;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_verQ;
      RD_fgdVA_verQ = object->Get<COMET::IRealDatum>("fgdVA_verQ");
      if( RD_fgdVA_verQ && i == fgdNum ){
   fgd->fgdVA_verQ = RD_fgdVA_verQ->GetValue();
      }
      else
        fgd->fgdVA_verQ = 0;
      
      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_verLayQ;
      RD_fgdVA_verLayQ = object->Get<COMET::IRealDatum>("fgdVA_verLayQ");
      if( RD_fgdVA_verLayQ && i == fgdNum ){
   fgd->fgdVA_verLayQ = RD_fgdVA_verLayQ->GetValue();
      }
      else
        fgd->fgdVA_verLayQ = 0;
      
      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_verNearQ;
      RD_fgdVA_verNearQ = object->Get<COMET::IRealDatum>("fgdVA_verNearQ");
      if( RD_fgdVA_verNearQ && i == fgdNum ){
   fgd->fgdVA_verNearQ = RD_fgdVA_verNearQ->GetValue();
      }
      else
        fgd->fgdVA_verNearQ = 0;
     
     COMET::IHandle<COMET::IRealDatum> RD_fgdVA_verNextNearQ;
      RD_fgdVA_verNextNearQ = object->Get<COMET::IRealDatum>("fgdVA_verNextNearQ");
      if( RD_fgdVA_verNextNearQ && i == fgdNum ){
   fgd->fgdVA_verNextNearQ = RD_fgdVA_verNextNearQ->GetValue();
      }
      else
        fgd->fgdVA_verNextNearQ = 0;
  
      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_verNextNextNearQ;
      RD_fgdVA_verNextNextNearQ = object->Get<COMET::IRealDatum>("fgdVA_verNextNextNearQ");
      if( RD_fgdVA_verNextNextNearQ && i == fgdNum ){
   fgd->fgdVA_verNextNextNearQ = RD_fgdVA_verNextNextNearQ->GetValue();
      }
      else
        fgd->fgdVA_verNextNextNearQ = 0;
  
      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_totalQ;
      RD_fgdVA_totalQ = object->Get<COMET::IRealDatum>("fgdVA_totalQ");
      if( RD_fgdVA_totalQ && i == fgdNum ){
   fgd->fgdVA_totalQ = RD_fgdVA_totalQ->GetValue();
      }
      else
        fgd->fgdVA_totalQ = 0;

      COMET::IHandle<COMET::IRealDatum> RD_fgdVA_totalCorrE;
      RD_fgdVA_totalCorrE = object->Get<COMET::IRealDatum>("fgdVA_totalCorrE");
      if( RD_fgdVA_totalCorrE && i == fgdNum ){
   fgd->fgdVA_totalCorrE = RD_fgdVA_totalCorrE->GetValue();
      }
      else
        fgd->fgdVA_totalCorrE = 0;
    }
  }      


}


//*************************************************************
void COMET::IGlobalReconModule::FillP0DInfo(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject  ){
//*************************************************************

  if (!object->UsesDetector(COMET::IReconBase::kP0D)) return;
  
  COMET::IHandle<COMET::IReconBase> p0dObject = ReconObjectUtils::GetConstituentInDetector(object, COMET::IReconBase::kP0D);  

  if (p0dObject && globalObject.NP0Ds < NMAXP0D){    
    
    IP0DObject* p0d = new((*(globalObject.P0D))[globalObject.NP0Ds++]) IP0DObject;
    
    // fill general information
    FillSubBaseObject(p0dObject, *p0d, 1);

    p0d->Cone          = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    p0d->Width         = DEFAULT_MAX;
    p0d->EDeposit      = DEFAULT_MAX;
      
    // The most  likely PID is  available through particle->GetParticle,
    // but it is also stored as one of the alternates.  These are sorted
    // in order of PID, so copy them directly into the TTree.
    COMET::IReconObjectContainer::const_iterator alternate;
    
    COMET::IHandle<COMET::IReconPID> pid = p0dObject;
    
    for(alternate = pid->GetAlternates().begin();
  alternate != pid->GetAlternates().end();
  alternate++){
      COMET::IHandle<COMET::IReconPID> alt = (*alternate);
      if(alt){
  p0d->ParticleId.push_back(alt->GetParticleId());
  p0d->PIDWeight.push_back(alt->GetPIDWeight());
      }
      else COMETWarn("IReconPID alternate is not a IReconPID - Ignored");
    }



    COMET::IHandle<COMET::IReconObjectContainer> conContainer = p0dObject->GetConstituents();
    if (conContainer){
      // The  PID variables are not saved in the main object but in its constituent
      COMET::IHandle<COMET::IReconBase> p0dObject2 = *(conContainer->begin());
      if (p0dObject2){ 

  // The cone angle is only accesible form the shower object
  p0d->Cone  = TrackingUtils::GetCone(p0dObject2);	
  p0d->EDeposit = TrackingUtils::GetEDeposit(p0dObject2);	

  if (GetObjectType(p0dObject).find("COMET::IReconShower") == std::string::npos) return; 

  // Compute the length and the width for showers
  COMET::IHandle<COMET::IHitSelection> hits = p0dObject2->GetHits();
  if (hits){
    
    double totalCharge = 0.0;
    double avgLen = 0.0;
    double avgLen2 = 0.0;
    
    for (COMET::IHitSelection::iterator h = hits->begin();h != hits->end();++h) {
      // Find the width.
      TVector3 diff = (*h)->GetPosition() - p0d->FrontPosition.Vect();
      TVector3 perp = diff - (diff*p0d->FrontDirection)*p0d->FrontDirection;
      if ((*h)->IsXHit()) {
        diff.SetY(0.0);
        perp.SetY(0.0);
      }
      else {
        diff.SetX(0.0);
        perp.SetX(0.0);
      }
      double q =  std::max(3.0,(*h)->GetCharge());
      p0d->Width += q*perp.Mag();
      double len = diff.Mag();
      avgLen += q*len;
      avgLen2 += q*len*len;
      totalCharge += q;
    }
    p0d->Width /= std::max(totalCharge,1.0);
    avgLen /= std::max(totalCharge,1.0);
    avgLen2 /= std::max(totalCharge,1.0);
    p0d->Length = avgLen2 - avgLen*avgLen;
    p0d->Length = std::sqrt(std::max(0.0,p0d->Length));
  }
      }
    }
  }


}

//*************************************************************
void COMET::IGlobalReconModule::FillSMRDInfo(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject  ){
//*************************************************************

  COMET::IReconBase::Status dets[4] ={COMET::IReconBase::kTopSMRD, COMET::IReconBase::kBottomSMRD, 
           COMET::IReconBase::kLeftSMRD, COMET::IReconBase::kRightSMRD};
  
  for (int i = 0;i<4;i++){
    if (!object->UsesDetector(dets[i])) continue;
    
    COMET::IHandle<COMET::IReconBase> smrdObject = ReconObjectUtils::GetConstituentInDetector(object, dets[i]);

  if (smrdObject && globalObject.NSMRDs < NMAXSMRD){    
      
      ISMRDObject* smrd = new((*(globalObject.SMRD))[globalObject.NSMRDs++]) ISMRDObject;

      // fill general information
      FillSubBaseObject(smrdObject, *smrd,i+1);
      
      //explicitly set SMRD EDeposit as a some of the charge of all hits
      smrd->EDeposit = DEFAULT_MAX;
      COMET::IHandle<COMET::IRealDatum> totalCharge =  smrdObject->Get< COMET::IRealDatum >("totalCharge");
      if (totalCharge)
  smrd->EDeposit = smrdObject->Get< COMET::IRealDatum >("totalCharge")->GetValue();
  
      // fill average  hit time
      smrd->avgtime = DEFAULT_MAX;
      COMET::IHandle<COMET::IRealDatum> smrdTimeRD = smrdObject->Get<COMET::IRealDatum>("averageHitTime");
      if (smrdTimeRD) 
  smrd->avgtime = smrdTimeRD->GetValue();
      
    }
  }      


}


//*************************************************************
void COMET::IGlobalReconModule::InitializeExtrapolationToDetectors(){
//*************************************************************


  std::vector<COMET::IGeomModuleBase*>::const_iterator it;

  fALLMODULES.resize(NDETSEXTRAP);

  // initialize all names
  for (int i=0;i<NDETSEXTRAP;i++)
    fALLMODULES[i]="";

  // Get TPC volumes
  std::vector<COMET::IGeomModuleBase*>  TPCMODULES = COMET::IGeomInfo::Get().TPC().GetModules();
  int index=0;
  for (it=TPCMODULES.begin();it!=TPCMODULES.end();it++){
    std::string name = (*it)->GetPath()+"/GasGap_0/Drift_0";
    fALLMODULES[index]=name;
    index++;
  }

  // Get FGD volumes 
  std::vector<COMET::IGeomModuleBase*>  FGDMODULES = COMET::IGeomInfo::Get().FGD().GetModules(); 
  index=3;
  for (it=FGDMODULES.begin();it!=FGDMODULES.end();it++){
    fALLMODULES[index]=(*it)->GetPath()+"/Active_0";
    index++;
  }

  // Get P0D volumes 
  std::vector<COMET::IGeomModuleBase*>  P0DMODULES = COMET::IGeomInfo::Get().P0D().GetModules();
  index=5;
  for (it=P0DMODULES.begin();it!=P0DMODULES.end();it++) {
    fALLMODULES[index]=(*it)->GetPath();
    index++;
  }
  
  // Get DsECAL volume
  const std::vector<COMET::IGeomModuleBase*>  ECALMODULES = COMET::IGeomInfo::Get().ECAL().GetModules();
  index=6;
  for (it=ECALMODULES.begin();it!=ECALMODULES.end();it++){
    if ((*it)->GetPath().find("DsECal") != std::string::npos){ 
      fALLMODULES[index]=(*it)->GetPath()+"/Active_0";
      index++;
    }
  }

  // Get SMRD volumes 
  std::string SMRD_name[6];
  std::string magnet_name = COMET::tman().GetSetup().GetVolumePath("Magnet_0");
  if (magnet_name != ""){
  
    SMRD_name[0] = magnet_name + "/RightClam_0/SMRD_0/Side_0";
    SMRD_name[1] = magnet_name + "/RightClam_0/SMRD_0/Top_0";
    SMRD_name[2] = magnet_name + "/RightClam_0/SMRD_0/Bottom_0";
    SMRD_name[3] = magnet_name + "/LeftClam_0/SMRD_0/Side_0";
    SMRD_name[4] = magnet_name + "/LeftClam_0/SMRD_0/Top_0";
    SMRD_name[5] = magnet_name + "/LeftClam_0/SMRD_0/Bottom_0";

    index=7;
    for (int i=0;i<6;i++){
      fALLMODULES[index]=SMRD_name[i];
      index++;
    }
  }

  // Get Barrel ECAL volumes
  index=13;
  for (it=ECALMODULES.begin();it!=ECALMODULES.end();it++){
    if ((*it)->GetPath().find("DsECal") == std::string::npos){ 
      fALLMODULES[index]=(*it)->GetPath()+"/Active_0";
      index++;
    }
  }


  //  const int NDETS = fALLMODULES.size();

  
  for (int i=0;i<NDETSEXTRAP;i++){
    std::string volname = fALLMODULES[i];
    if (volname=="") continue;
    const Volume& vol = COMET::gman().GetSetup().volume(volname);
    const surface_vector& surfs = vol.surfaces("outer");
    stc_tools::extend_vector<Surface*>(surfs,fDetSurfaces[i]);    

    fDetIndex[vol.name("setup_name")] = i;
  }

  if (debug_extrap){

    std::cout << "list of volumes: " << std::endl;
    for (int i=0;i<NDETSEXTRAP;i++){
      std::cout << i << " --> " << fALLMODULES[i] << std::endl;
    }

    std::cout << "fDetIndex" << std::endl;
    //  std::map<std::string, int>::const_iterator it2;  
    //  dict::dictionary<int>::const_iterator it2;  
    for (unsigned int it2=0;it2<fDetIndex.keys().size();it2++){
      std::cout << fDetIndex.keys()[it2] << " " << fDetIndex[fDetIndex.keys()[it2]] << std::endl; 
      int i = fDetIndex[fDetIndex.keys()[it2]];
      for (unsigned int j=0;j<fDetSurfaces[i].size();j++)
  std::cout << "    " << j << ":  " << (fDetSurfaces[i])[j]->name("setup_name") << std::endl;
    }
  }
}

//*************************************************************
void COMET::IGlobalReconModule::FillExtrapolationToDetectors(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject  ){
//*************************************************************

  if (debug_extrap){
    std::cout << "--------------------------------------------------" << std::endl;
    std::cout << "IGlobalReconModule:: Extrapolate to subdetectors: " << std::endl;
    std::cout << "--------------------------------------------------" << std::endl;
  }


  // Initialize variables
  for (int i=0;i<NDETSEXTRAP;i++){
    fPassedDetector[i]=false;
    globalObject.EntranceOK[i]=0;
    globalObject.ExitOK[i]=0;

    globalObject.EntrancePosition[i]=TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    globalObject.ExitPosition[i]=TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);

    globalObject.EntranceDirection[i]=TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    globalObject.ExitDirection[i]=TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);

    globalObject.EntranceMomentum[i]=DEFAULT_MAX;
    globalObject.ExitMomentum[i]=DEFAULT_MAX;

    globalObject.EntrancePositionErr[i]=TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    globalObject.ExitPositionErr[i]=TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);

    globalObject.EntranceDirectionErr[i]=TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    globalObject.ExitDirectionErr[i]=TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);

    globalObject.EntranceMomentumErr[i]=DEFAULT_MAX;
    globalObject.ExitMomentumErr[i]=DEFAULT_MAX;
  }

  // convert the track to a RecPack trajectory
  Trajectory traj;
  bool ok = COMET::converter().TReconBase_to_Trajectory(object, traj);
  if (!ok) return;

  if (!traj.status("fitted")){
    return;
  }

  if (traj.state(traj.first_fitted_node()).vector()[5] < 0){
    traj.sort_nodes(RP::z,-1);
    if (debug_extrap)
      std::cout << "Negative sense: sort nodes in reverse z order" << std::endl;
  }

  if (debug_extrap)
    std::cout << "Trajectory: " << traj << std::endl;



  // extrapolate forwards
  FillExtrapolationToDetectors(object, globalObject, traj, 1);  

  // extrapolate backwards
  FillExtrapolationToDetectors(object, globalObject, traj, -1);  

  if (debug_extrap){
    std::cout << "---------------------------------------------------------" << std::endl;
    std::cout << "IGlobalReconModule:: Extrapolate to subdetectors results:" << std::endl;
    std::cout << "---------------------------------------------------------" << std::endl;

    std::cout << object << std::endl;
  
    //    for (unsigned int i =0; i<globalObject.EntrancePosition.size(); i++){
    for (int i =0; i<NDETSEXTRAP; i++){
      if (globalObject.EntranceOK[i])
  std::cout << i 
      << "   x = " << globalObject.EntrancePosition[i].X()
      << "\t y = " << globalObject.EntrancePosition[i].Y()
      << "\t z = " << globalObject.EntrancePosition[i].Z()
      << "\t ex = " << globalObject.EntrancePositionErr[i].X()
      << "\t ey = " << globalObject.EntrancePositionErr[i].Y()
      << "\t ez = " << globalObject.EntrancePositionErr[i].Z()
      << std::endl;
      if (globalObject.ExitOK[i])
  std::cout << i 
      << "   x = " << globalObject.ExitPosition[i].X() 
      << "\t y = " << globalObject.ExitPosition[i].Y()
      << "\t z = " << globalObject.ExitPosition[i].Z()
      << "\t ex = " << globalObject.ExitPositionErr[i].X()
      << "\t ey = " << globalObject.ExitPositionErr[i].Y()
      << "\t ez = " << globalObject.ExitPositionErr[i].Z()
      << std::endl;
      
      if (globalObject.EntranceOK[i])
  std::cout << i 
      << "   ux = " << globalObject.EntranceDirection[i].X() 
      << "\t uy = " << globalObject.EntranceDirection[i].Y()
      << "\t uz = " << globalObject.EntranceDirection[i].Z()
      << "\t ex = " << globalObject.EntranceDirectionErr[i].X()
      << "\t ey = " << globalObject.EntranceDirectionErr[i].Y()
      << "\t ez = " << globalObject.EntranceDirectionErr[i].Z()
      << std::endl;
      if (globalObject.ExitOK[i])
  std::cout << i 
      << "   ux = " << globalObject.ExitDirection[i].X() 
      << "\t uy = " << globalObject.ExitDirection[i].Y()
      << "\t uz = " << globalObject.ExitDirection[i].Z()
      << "\t ex = " << globalObject.ExitDirectionErr[i].X()
      << "\t ey = " << globalObject.ExitDirectionErr[i].Y()
      << "\t ez = " << globalObject.ExitDirectionErr[i].Z()
      << std::endl;
      
      if (globalObject.EntranceOK[i])
  std::cout << i 
      << "   P0 = " << (globalObject.EntranceMomentum[i])
      << "   eP0 = " << (globalObject.EntranceMomentumErr[i])
      << std::endl;
      if (globalObject.ExitOK[i])
  std::cout << i 
      << "   P1 = " << (globalObject.ExitMomentum[i]) 
      << "   eP1 = " << (globalObject.ExitMomentumErr[i]) 
      << std::endl;
      
      
    }
    
    
  }


}

//*************************************************************
void COMET::IGlobalReconModule::FillExtrapolationToDetectors(COMET::IHandle<COMET::IReconBase> object, IGlobalPID& globalObject, Trajectory& traj, int sense  ){
//*************************************************************

  // 1. Take the first state in the trajectory (taking into account the sense)
  int ifirst = traj.first_fitted_node();
  State state  = traj.state(ifirst);

  if (debug_extrap){
    std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
    std::cout << "1. First state: " << state << std::endl;

  }

  // 2. Find the volume associated to the first state

  const Volume* volume = NULL;
  EVector position = fEquation.position(state,0.);
  bool ok = COMET::gman().GetSetup().volume(position,volume);
  if (!ok){
    if (debug_extrap)
      std::cout << "Volume not found. Position:  " << position << std::endl;
    return;
  }

  int current_det = -1;
  for (unsigned int i=0;i<fALLMODULES.size();i++){
    if (volume->name("setup_name").find(fALLMODULES[i]) != std::string::npos)
      current_det = i;
  }

  if (current_det<0){
    if (debug_extrap){
      std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
      std::cout << "2. Current volume not in list, name = " << volume->name("setup_name") << std::endl;
    }
    return;
  }

  if (debug_extrap){
    std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
    std::cout << "2. Current volume: index = " << current_det << ", name = " << volume->name("setup_name") << std::endl;
  }

  fPassedDetector[current_det] = false;
  int iexit=-1;
  bool in_found = true;

  while (volume){
    
    if (debug_extrap){
      std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
      std::cout << "3. Exit from volume: " << volume->name("setup_name")<< std::endl;     
      std::cout << "   ifirst: " << ifirst << ", iexit = " << iexit << ", in found = " << in_found <<  std::endl; 
    }    

    // 3. Loop over states until we exit from the volume. Only when there is an state inside the current volume (in_found=true)
    if (in_found){
      // by default the exit state is the first one
      iexit = ifirst;
      if (sense==1){
  int iprev=ifirst;
  for (unsigned int i=ifirst+1;i<traj.nodes().size();i++){
    if (traj.node(i).status("fitted")){
      position = fEquation.position( traj.state(i),0.);
      bool ok = volume->is_inside(position);
      if (!ok){
        // this is the first state outside the volume. Take the previous as exit state
        iexit = iprev;
        break;
      }
      iprev=i;
    }
  }
      }
      if (debug_extrap)
  std::cout << "   Last state in volume: " << iexit << std::endl;

      // take the state
      state = traj.state(iexit);

      if (debug_extrap)
  std::cout << state << std::endl;
    }
    else{
      if (debug_extrap)
  std::cout << "   No state inside volume. Propagate to exit surface from the entrance state " << std::endl;      
    }

    // 4. Find the closest surface to the last state in the volume (or the previous extrapolation if there is no state in the volume)
    surface_vector exit_surfaces;
    stc_tools::extend_vector<Surface*>(fDetSurfaces[current_det],exit_surfaces);
    
    double end_length = 1000000*sense;
    double delta_length;
    const Surface* exit_surf;
    
    ok = COMET::rpman("globalAnalysis").navigation_svc().next_surface(state,end_length, exit_surfaces, delta_length, exit_surf);

    if (ok){
      if (debug_extrap){
  std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
  std::cout << "4. Exit surface. Length: " << delta_length << std::endl;
  std::cout << "                 name:   " << exit_surf->name("setup_name") << std::endl;
  std::cout << "                 volume: " << volume->name("setup_name") << std::endl;
      }
    }
    else{
      if (debug_extrap){
  std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
  std::cout << "4. Exit surface. Fail in finding exit surface. stop !!! " << std::endl;
      }
      return;
    }

    // 5. propagate to the surface
    ok = COMET::rpman("globalAnalysis").navigation_svc().propagate(0.1*sense,state);
    ok = COMET::rpman("globalAnalysis").navigation_svc().propagate(*exit_surf, state, delta_length);

    if (ok){
      // take the position as the exit position
      EVector exitPosition  = fEquation.position(state,0.);
      EVector exitDirection = fEquation.direction(state,0.);
      double exitMomentum   = DEFAULT_MAX;
      if (state.vector()[6] !=0) 
  exitMomentum = fabs(1/state.vector()[6]);


      TLorentzVector posErr(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
      if (state.matrix()[0][0]>0 && state.matrix()[1][1]>0 && state.matrix()[2][2]>0)
  posErr = TLorentzVector(sqrt(state.matrix()[0][0]), sqrt(state.matrix()[1][1]), sqrt(state.matrix()[2][2]),0) ;

      TVector3 dirErr(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX); 
      if (state.matrix()[3][3]>0 && state.matrix()[4][4]>0 && state.matrix()[5][5]>0)
  dirErr = TVector3(sqrt(state.matrix()[3][3]), sqrt(state.matrix()[4][4]), sqrt(state.matrix()[5][5])); 

      double exitMomentumErr   = DEFAULT_MAX;
      if (state.matrix()[6][6] >0 && state.vector()[6] !=0) 
  exitMomentumErr = fabs(1/pow(state.vector()[6],2))*sqrt(state.matrix()[6][6]);

      
      if (debug_extrap){
  std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
  std::cout << "5. Exit position:  " << print(exitPosition) << std::endl;
  std::cout << "   Exit direction: " << print(exitDirection) << std::endl;
  std::cout << "   Exit momentum:  " << exitMomentum << std::endl;
      }            
     
      TLorentzVector pos;
      ok = COMET::converter().EVector3_to_TVector(exitPosition, pos);

      TLorentzVector dir4;
      ok = COMET::converter().EVector3_to_TVector(exitDirection, dir4);
      TVector3 dir = dir4.Vect();


      if (sense==1){
  globalObject.ExitOK[current_det] = 1;
  globalObject.ExitPosition[current_det] = pos;
  globalObject.ExitDirection[current_det] = dir;
  globalObject.ExitMomentum[current_det] = exitMomentum;
  globalObject.ExitPositionErr[current_det] = posErr;
  globalObject.ExitDirectionErr[current_det] = dirErr;
  globalObject.ExitMomentumErr[current_det] = exitMomentumErr;
      }
      else{
  globalObject.EntranceOK[current_det] = 1;
  globalObject.EntrancePosition[current_det] = pos;
  globalObject.EntranceDirection[current_det] = dir;
  globalObject.EntranceMomentum[current_det] = exitMomentum;
  globalObject.EntrancePositionErr[current_det] = posErr;
  globalObject.EntranceDirectionErr[current_det] = dirErr;
  globalObject.EntranceMomentumErr[current_det] = exitMomentumErr;

      }


    }
    else{
      if (debug_extrap){
  std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
  std::cout << "5. Exit position: " << "surface not intersected !!!" << std::endl;
      }
    }

    fPassedDetector[current_det] = true;
    
    // 6. Find the closest surface to the exit state among all other volumes

    // list of possible entrance surfaces 
    surface_vector entrance_surfaces;
    for (int det=0;det<NDETSEXTRAP;det++){
      if (!fPassedDetector[det]){
  stc_tools::extend_vector<Surface*>(fDetSurfaces[det],entrance_surfaces);    
      }
    }
    
    const Surface* entrance_surf;
    ok = COMET::rpman("globalAnalysis").navigation_svc().next_surface(state,end_length, entrance_surfaces, delta_length, entrance_surf);
    
    if (!ok){
      if (debug_extrap)
  std::cout << "No suface intersected " << std::endl;
      break;
    }

    // find the parent volume 
    const Volume* volume_try;
    ok = COMET::gman().GetSetup().volume(*entrance_surf,volume_try);
    
    if (debug_extrap){
      std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
      std::cout << "6. Entrance surface. Length: " << delta_length << std::endl;
      std::cout << "                     name:   " << entrance_surf->name("setup_name") << std::endl;
      std::cout << "                     volume: " << volume_try->name("setup_name") << std::endl;
    }
    
    //7.  Check if the next state in the trajectory is in that volume    
    in_found = false;
    if (sense==1){
      if ((unsigned int)(iexit+1) < traj.nodes().size()){
  for (unsigned int i=iexit+1;i<traj.nodes().size();i++){
    if (traj.node(i).status("fitted")){
      position = fEquation.position( traj.state(i),0.);
      bool ok = volume_try->is_inside(position);
      if (ok){
        ifirst = i;
        in_found = true;
      }
      break;	
    }
  }
      }
    }
    
    if (debug_extrap){
      std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
      std::cout << "7. Enter in volume: " << volume_try->name("setup_name") << std::endl;
      std::cout << "   ifirst: " << ifirst << ", iexit = " << iexit << ", in found = " << in_found <<  std::endl; 
    }

    //    iexit=ifirst
   
    if (debug_extrap){
      if (in_found)
  std::cout << "   Next state inside volume: " << ifirst << std::endl;
      else
  std::cout << "   No state inside volume. Propagate to entrance surface from the exit of the previous volume " << std::endl;
    }

    // 8.  propagate to the surface

    if (in_found){
      // from the first state in the current volume
      state = traj.state(ifirst);
    }

    ok = COMET::rpman("globalAnalysis").navigation_svc().propagate(0.1*sense,state);
    ok = COMET::rpman("globalAnalysis").navigation_svc().propagate(*entrance_surf, state, delta_length);

    if (ok){
      volume = volume_try;
      current_det = -1;
      for (unsigned int i=0;i<fALLMODULES.size();i++){
  if (volume->name("setup_name").find(fALLMODULES[i]) != std::string::npos)
    current_det = i;
      }

      // take the position as the entrance position
      EVector entrancePosition  = fEquation.position(state,0.);
      EVector entranceDirection = fEquation.direction(state,0.);
      double entranceMomentum   = DEFAULT_MAX;

      if (state.vector()[6] !=0) 
  entranceMomentum   = fabs(1/state.vector()[6]);

      TLorentzVector posErr(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
      if (state.matrix()[0][0]>0 && state.matrix()[1][1]>0 && state.matrix()[2][2]>0)
  posErr = TLorentzVector(sqrt(state.matrix()[0][0]), sqrt(state.matrix()[1][1]), sqrt(state.matrix()[2][2]),0) ;

      TVector3 dirErr(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX); 
      if (state.matrix()[3][3]>0 && state.matrix()[4][4]>0 && state.matrix()[5][5]>0)
  dirErr = TVector3(sqrt(state.matrix()[3][3]), sqrt(state.matrix()[4][4]), sqrt(state.matrix()[5][5])); 

      double entranceMomentumErr   = DEFAULT_MAX;
      if (state.matrix()[6][6] >0 && state.vector()[6] !=0) 
  entranceMomentumErr = fabs(1/pow(state.vector()[6],2))*sqrt(state.matrix()[6][6]);

      
   
      if (debug_extrap){   
  std::cout << "8. Entrance position:  " << print(entrancePosition) << std::endl;
  std::cout << "   Entrance direction: " << print(entranceDirection) << std::endl;
  std::cout << "   Entrance momentum:  " << entranceMomentum << std::endl;
      }            

      TLorentzVector pos;
      ok = COMET::converter().EVector3_to_TVector(entrancePosition, pos);

      TLorentzVector dir4;
      ok = COMET::converter().EVector3_to_TVector(entranceDirection, dir4);
      TVector3 dir = dir4.Vect();

      if (sense==1){
  globalObject.EntranceOK[current_det] = 1;
  globalObject.EntrancePosition[current_det] = pos;
  globalObject.EntranceDirection[current_det] = dir;
  globalObject.EntranceMomentum[current_det] = entranceMomentum;
  globalObject.EntrancePositionErr[current_det] = posErr;
  globalObject.EntranceDirectionErr[current_det] = dirErr;
  globalObject.EntranceMomentumErr[current_det] = entranceMomentumErr;
      }
      else{
  globalObject.ExitOK[current_det] = 1;
  globalObject.ExitPosition[current_det] = pos;
  globalObject.ExitDirection[current_det] = dir;
  globalObject.ExitMomentum[current_det] = entranceMomentum;
  globalObject.ExitPositionErr[current_det] = posErr;
  globalObject.ExitDirectionErr[current_det] = dirErr;
  globalObject.ExitMomentumErr[current_det] = entranceMomentumErr;
      }
    }
    else{
      if (debug_extrap){
  std::cout << "Extrapolate to subdetectors: " <<  std::endl;   
  std::cout << "8. Entrance position: " << "surface not intersected !!!" << std::endl;
      }
    }

    // 9. Repit 3 to 8 for all other volumes. 

  }


}


//***********************************************************************************
bool COMET::IGlobalReconModule::FillVertex(COMET::IHandle<COMET::IReconVertex> vertex, bool primary){
//***********************************************************************************

  // Fill a reconstructed vertex and associated tracks
    
  //Maximum number of secondary vertices in spill
  if(!primary && fNSVertices == NSVERTICES) return false;
 
  COMET::IHandle<COMET::IVertexState> vs = vertex->GetState();
  if(!vs) return false;
 
  IGlobalVertex* gVertex = new((*fVertices)[fNVertices++]) IGlobalVertex;

  gVertex->NTrueVertices = 0;

  //Primary
  if(primary){
    fPVInd = fNVertices - 1; //index of last added primary
    gVertex->PrimaryIndex = -1; //
  }
  else{
    //For secondaries set the index of the corresponding primary vertex 
    gVertex->PrimaryIndex = fPVInd; 
    fNSVertices ++;
  }
  gVertex->Position      = vs->GetPosition();
  gVertex->Variance      = vs->GetPositionVariance();  
  gVertex->AlgorithmName = vertex->GetAlgorithmName();
  gVertex->Status        = vertex->GetStatus();
  gVertex->Quality       = vertex->GetQuality();
  gVertex->NDOF          = vertex->GetNDOF();

  COMET::IHandle<COMET::IReconObjectContainer> tracksAssociated = vertex->GetConstituents();
  if(tracksAssociated->size() == 0) return true;

  std::map<COMET::IHandle<COMET::IReconBase>, COMET::IHandle<COMET::IReconObjectContainer> >::iterator cit;

  //Loop over associated tracks
  for(COMET::IReconObjectContainer::iterator assocIter = tracksAssociated->begin();
      assocIter != tracksAssociated->end(); assocIter++){

    COMET::IHandle<COMET::IReconBase> track = *assocIter;
    if(!track) continue;

    COMET::IHandle<COMET::IReconState> state = track->GetState();
    if(!state) continue;

    //The original track stored as constituent for Kalman vertices
    COMET::IHandle<COMET::IReconBase> otrack;
    if(gVertex->AlgorithmName == "Kalman"){
      COMET::IHandle<COMET::IReconObjectContainer> trackC = track->GetConstituents();
      if(trackC->size() == 0) continue;
      otrack = trackC->front();
      if(!otrack) continue;
    }

    double mom = TrackingUtils::GetMomentum(state);
    TVector3 mdir = TrackingUtils::GetDirection(state);
    mdir *= mom;
    int charge = (int)TrackingUtils::GetCharge(track);
    double quality = track->GetQuality();

    //Fill track info
    if(gVertex->NConstituents == NCONSTITUENTS) break;
    IVertexConstituent* vConst = new((*(gVertex->Constituents))[gVertex->NConstituents++]) IVertexConstituent;

    vConst->Charge = charge;
    vConst->Quality = quality;
    vConst->Momentum = mdir;

    //Get PID info
    int PID = -1;
    bool inMap = false;
    std::map<COMET::IHandle<COMET::IReconBase>, int >::iterator it, miter;
    //Use original track for Kalman vertices
    TLorentzVector pos = TrackingUtils::GetPosition(state);
    if(otrack){
      COMET::IHandle<COMET::IReconState> ostate = otrack->GetState();
      mom = TrackingUtils::GetMomentum(ostate);
      pos = TrackingUtils::GetPosition(ostate);
    }
    //Global vertex
    for(it = fGlobalIndexMap.begin(); it != fGlobalIndexMap.end(); it ++){
      COMET::IHandle<COMET::IReconState> gstate = (it->first)->GetState();
      double gmom = TrackingUtils::GetMomentum(gstate);
      TLorentzVector gpos = TrackingUtils::GetPosition(gstate);
      if(gpos == pos && gmom == mom){
  miter = fGlobalIndexMap.find(it->first);
  if(miter != fGlobalIndexMap.end()){
    inMap = true;
    break;
  }
      }
    }
    //Tracker vertex
    if(!inMap){
      if(otrack) miter = fTrackerGlobalIndexMap.find(otrack);
      else miter = fTrackerGlobalIndexMap.find(track);
      if(miter != fTrackerGlobalIndexMap.end()) inMap = true;
    }
    //Set PID
    if(inMap) PID = miter->second;
    vConst->PID = PID;
  }

return true;

}


//*************************************************************
void COMET::IGlobalReconModule::MatchTrueVertex(COMET::IHandle<COMET::IG4PrimaryVertexContainer> g4PrimVert){
//*************************************************************

  IGlobalVertex* gVertex = (IGlobalVertex*)(*fVertices)[fNVertices-1];

  std::vector<TLorentzVector> vertexPos;
  std::vector< std::vector<int> > vertexPIDs;

  int nConst = gVertex->NConstituents;

  //Fill true vertices
  //Loop over constituents
  for(int i = 0; i < nConst; i ++){

    if(gVertex->NTrueVertices == NCONSTITUENTS) break;
    IVertexConstituent* vConst = (IVertexConstituent*)(*(gVertex->Constituents))[i];
    int this_PID = vConst->PID;

    IGlobalPID* globalObject = (IGlobalPID*)(*fPIDs)[this_PID];
    ITrueParticle globalTrue = globalObject->TrueParticle;

    TLorentzVector vPos = globalTrue.Vertex.Position;
    bool inList = false;
    for(int j = 0; j < (int)vertexPos.size(); j ++){
      if(vertexPos[j] == vPos) inList = true;
    }
    if(inList) continue;
    vertexPos.push_back(vPos);

    //Loop over all PIDs, and store those that come from this vertex
    std::vector<int> IDs;
    for(int k = 0; k < fNPIDs; k ++){
      IGlobalPID* thisGlobalObject = (IGlobalPID*)(*(fPIDs))[k];
      ITrueParticle thisGlobalTrue = thisGlobalObject->TrueParticle;
      TLorentzVector thisVPos = thisGlobalTrue.Vertex.Position;
      if((thisVPos - vPos).Mag() < 0.01) IDs.push_back(k);
    }
    vertexPIDs.push_back(IDs);

    ITrueVertex* tVertex = new((*(gVertex->TrueVertices))[gVertex->NTrueVertices++]) ITrueVertex;
    *tVertex = globalTrue.Vertex;
  }

  //Set efficiencies and purities
  //Loop over the found true vertices
  for(int i = 0; i < gVertex->NTrueVertices; i ++){

    std::vector<int> IDs = vertexPIDs[i];
    int nIDs = (int)IDs.size();
    int nMatch = 0;

    //Loop over vertex constituents
    for(int j = 0; j < nConst; j ++){
      IVertexConstituent* vConst = (IVertexConstituent*)(*(gVertex->Constituents))[j];
      int this_PID = vConst->PID;
      //Check if they match true vertex PID list
      for(int k = 0; k < nIDs; k ++){
  if(this_PID == IDs[k]) nMatch ++;
      }
    }

    ITrueVertex* tVertex = (ITrueVertex*)(*gVertex->TrueVertices)[i];

    tVertex->Pur = (double)nMatch/nConst;
    tVertex->Eff = (double)nMatch/nIDs;
  }
}


//*************************************************************
void COMET::IGlobalReconModule::FillTrueParticle( COMET::IHandle<COMET::IG4Trajectory> G4track, double pur, double eff, COMET::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, DEFAULT_MAX, DEFAULT_MAX, trueParticle.Vertex);

  }
  else{

    trueParticle.ID    = -1;
   
    trueParticle.Pur = DEFAULT_MAX;
    trueParticle.Eff = DEFAULT_MAX;
        
    COMET::IG4PrimaryVertex G4vertex;
    FillTrueVertex(false, G4vertex, DEFAULT_MAX, DEFAULT_MAX, trueParticle.Vertex);
  }

}

//*************************************************************
COMET::IHandle<COMET::IG4Trajectory> COMET::IGlobalReconModule::GetParent(COMET::IHandle<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>();

  /*
  // create a map with trajectories ordered in descending kinetic energy
  if(trajectories){
    COMET::IG4TrajectoryContainer::iterator trajIter;
    for(trajIter = trajectories->begin();trajIter != trajectories->end(); ++trajIter) {
      if ((*trajIter)->GetTrackId() == G4track->GetParentId()) return *trajIter;
    }
  }

  return COMET::IHandle<COMET::IG4Trajectory>();
  */
}

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

  if (found){
    vertex.Position = G4vertex.GetPosition();
    vertex.ID = G4vertex.GetInteractionNumber();
    vertex.Pur = pur;
    vertex.Eff = eff;
    }
  else{
    vertex.Position = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    vertex.ID = -1;
    vertex.Pur = DEFAULT_MAX;
    vertex.Eff = DEFAULT_MAX;
  }

}

//*************************************************************
COMET::IHandle<COMET::IG4PrimaryVertex> COMET::IGlobalReconModule::GetG4Vertex(const COMET::IReconBase&, double& pur, double& eff){
//*************************************************************

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

}


//*************************************************************
bool  COMET::IGlobalReconModule::GetG4Vertex(COMET::IHandle<COMET::IG4Trajectory> G4track, COMET::IG4PrimaryVertex&  G4vertex){
//*************************************************************

  // 1. Find the primary particle associated to this G4 trajectory

  const COMET::ICOMETEvent& event = *(COMET::IEventFolder::GetCurrentEvent());  
  COMET::IHandle<COMET::IG4TrajectoryContainer> trajectories = event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");
  if (!trajectories) return false;
  
  int nextparent;
  int count = 0;
  int primaryID = G4track->GetTrackId();

  while((primaryID > 0) && ((nextparent = (*trajectories)[primaryID].GetParentId()) > 0) )  {
    primaryID = nextparent;
    if(count++ > 100000000)     {
      // just to make sure it doesn't go into an infinite loop somehow                                                                                                  
      COMETError("IGlobalReconModule: Particle has more than 100000000 levels of parents!");
      throw COMET::EoaAnalysisInfiniteLoop2();
    }
  }

  // 2. Find the primary vertex associated to the primary particle
  COMET::IHandle<COMET::IG4PrimaryVertexContainer>   vertices  = event.Get<COMET::IG4PrimaryVertexContainer>("truth/G4PrimVertex00");
  if (!vertices) return false;

  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() == primaryID){
  G4vertex = *it1;
  return true;
      }
    }
  }
  
  return false;

}


//*************************************************************
bool  COMET::IGlobalReconModule::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::IGlobalReconModule::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++){
      nHits += GetNumberOfHits(*it);
    }
  }

  return nHits;

}

//*****************************************************
int COMET::IGlobalReconModule::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::kP0D)    ) det=det*10+6;
  if (object->UsesDetector(COMET::IReconBase::kDSECal) ) det=det*10+7;
  if (object->UsesDetector(COMET::IReconBase::kTopSMRD)   ) det=det*10+8;
  if (object->UsesDetector(COMET::IReconBase::kBottomSMRD)) det=det*10+8;
  if (object->UsesDetector(COMET::IReconBase::kLeftSMRD)  ) det=det*10+8;
  if (object->UsesDetector(COMET::IReconBase::kRightSMRD) ) det=det*10+8;
  if (object->UsesDetector(COMET::IReconBase::kTopPECal)    ) det=det*10+9;   
  if (object->UsesDetector(COMET::IReconBase::kBottomPECal) ) det=det*10+9;   
  if (object->UsesDetector(COMET::IReconBase::kLeftPECal)   ) det=det*10+9;   
  if (object->UsesDetector(COMET::IReconBase::kRightPECal)  ) det=det*10+9;   
  if (object->UsesDetector(COMET::IReconBase::kTopTECal)    ) det=det*10+9;   
  if (object->UsesDetector(COMET::IReconBase::kBottomTECal) ) det=det*10+9;   
  if (object->UsesDetector(COMET::IReconBase::kLeftTECal)   ) det=det*10+9;   
  if (object->UsesDetector(COMET::IReconBase::kRightTECal)  ) det=det*10+9;   


  return det;

}


//*****************************************************
void COMET::IGlobalReconModule::FillDetectorUsed(COMET::IHandle<COMET::IReconBase> object, bool dets[]){
//*****************************************************

  for (int i=0;i<NDETSUSED;i++)
    dets[i]=false;

  if (object->UsesDetector(COMET::IReconBase::kTPC1)   ) dets[0]=true;
  if (object->UsesDetector(COMET::IReconBase::kTPC2)   ) dets[1]=true;
  if (object->UsesDetector(COMET::IReconBase::kTPC3)   ) dets[2]=true;
  if (object->UsesDetector(COMET::IReconBase::kFGD1)   ) dets[3]=true;
  if (object->UsesDetector(COMET::IReconBase::kFGD2)   ) dets[4]=true;
  if (object->UsesDetector(COMET::IReconBase::kP0D)    ) dets[5]=true;
  if (object->UsesDetector(COMET::IReconBase::kDSECal) ) dets[6]=true;
  if (object->UsesDetector(COMET::IReconBase::kTopSMRD)   ) dets[7]=true;
  if (object->UsesDetector(COMET::IReconBase::kBottomSMRD)) dets[8]=true;
  if (object->UsesDetector(COMET::IReconBase::kLeftSMRD)  ) dets[9]=true;
  if (object->UsesDetector(COMET::IReconBase::kRightSMRD) ) dets[10]=true;
  if (object->UsesDetector(COMET::IReconBase::kTopPECal)  )    dets[11]=true;
  if (object->UsesDetector(COMET::IReconBase::kBottomPECal)  ) dets[12]=true;
  if (object->UsesDetector(COMET::IReconBase::kLeftPECal)  )   dets[13]=true;
  if (object->UsesDetector(COMET::IReconBase::kRightPECal)  )  dets[14]=true;
  if (object->UsesDetector(COMET::IReconBase::kTopTECal)  )    dets[15]=true;
  if (object->UsesDetector(COMET::IReconBase::kBottomTECal)  ) dets[16]=true;
  if (object->UsesDetector(COMET::IReconBase::kLeftTECal)  )   dets[17]=true;
  if (object->UsesDetector(COMET::IReconBase::kRightTECal)  )  dets[18]=true;


}

//*****************************************************
void COMET::IGlobalReconModule::FillTPCOther(COMET::ICOMETEvent& event){
//*****************************************************


  COMET::IHandle<COMET::IAlgorithmResult> ReconTPC = event.GetFit("ReconTPC");
  if (ReconTPC){ 
    COMET::IHandle<COMET::IReconObjectContainer> otherObjects = ReconTPC->GetResultsContainer("other");
    if (otherObjects){  
      
      for(COMET::IReconObjectContainer::iterator tt = otherObjects->begin();
    tt != otherObjects->end(); tt++) {
  
  COMET::IHandle<COMET::IReconBase> other = *tt;

  if ( other &&  fNTPCOthers < NMAXTPCOTHER){    
    ITPCOtherObject* tpcOther = new((*fTPCOthers)[fNTPCOthers++]) ITPCOtherObject;



    tpcOther->Detector  = GetDetectorNumber(other);
    tpcOther->Chi2      = other->GetQuality();
    tpcOther->EDeposit = TrackingUtils::GetEDeposit(other);
    tpcOther->Charge   = TrackingUtils::GetCharge(other);	
    tpcOther->NHits    = DEFAULT_MAX;
    if (other->GetHits())
      tpcOther->NHits    = other->GetHits()->size();
    
    // Set the back and front position, variance and direction
    COMET::IHandle<COMET::IReconState> firstState = TrackingUtils::GetFirstState(*other);
    if (firstState){
      tpcOther->FrontPosition       = TrackingUtils::GetPosition(firstState);
      tpcOther->FrontDirection      = TrackingUtils::GetDirection(firstState);     
      tpcOther->Momentum       = TrackingUtils::GetMomentum(firstState);
    }
    else{
      tpcOther->FrontPosition       = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
      tpcOther->FrontDirection      = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
      tpcOther->Momentum       = DEFAULT_MAX; 
    }
    
    COMET::IHandle<COMET::IReconState> lastState  = TrackingUtils::GetLastState(*other);
    if (lastState){
      tpcOther->BackPosition       = TrackingUtils::GetPosition(lastState);
      tpcOther->BackDirection      = TrackingUtils::GetDirection(lastState);     
    }
    else{
      tpcOther->BackPosition       = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
      tpcOther->BackDirection      = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    }
    
    // Fill the true particle when it is found
    double pur, eff;
    COMET::IHandle<COMET::IG4Trajectory> G4track = GetG4Trajectory(*other,pur,eff);
    FillTrueParticle(G4track, pur, eff, tpcOther->TrueParticle);
    
  }
      }
    }
  }


}


//*****************************************************
void COMET::IGlobalReconModule::FillUnusedHits(COMET::ICOMETEvent& event){
//*****************************************************
  

  // Now add TPC hits
  COMET::IHandle<COMET::IAlgorithmResult> ReconTPC = event.GetFit("ReconTPC");
  if (ReconTPC){ 
    // Get TPC unused hits
    COMET::IHandle<COMET::IHitSelection> tpcHits = ReconTPC->GetHitSelection("unused");
    if (tpcHits){  
      fNTPCUnused = tpcHits->size();
    }
  }


  // Now add P0D hits
  COMET::IHandle<COMET::IAlgorithmResult> ReconP0D = event.GetFit("ReconP0D");
  if (ReconP0D){ 
    // Get P0D unused hits
    COMET::IHandle<COMET::IHitSelection> p0dHits = ReconP0D->GetHitSelection("unused");
    if (p0dHits){  
      for(COMET::IHitSelection::iterator hit = p0dHits->begin(); hit != p0dHits->end(); hit++) {
  IGlobalHit* gHit = new((*fP0DUnused)[fNP0DUnused++]) IGlobalHit;      
  FillGlobalHit(*hit, *gHit);
      }
    }
  }

  COMET::IHandle<COMET::IAlgorithmResult> ReconSMRD = event.GetFit("ReconSMRD");
  if (ReconSMRD){ 
    // Get SMRD unused hits
    COMET::IHandle<COMET::IHitSelection> smrdHits = ReconSMRD->GetHitSelection("unused");
    if (smrdHits){  
      for(COMET::IHitSelection::iterator h = smrdHits->begin(); h != smrdHits->end(); h++) {
        
  //fill info for smrd unused hit
  ISMRDHit* smrdHit=NULL;
        smrdHit = new((*fSMRDUnused)[fNSMRDUnused++]) ISMRDHit;      
  FillSmrdHit(*h, *smrdHit);  
  
  if (COMET::IGeomInfo::SMRD().IsSMRDTopWall((*h))) {
    fNSMRDTopUnused++;
  } 
  else if (COMET::IGeomInfo::SMRD().IsSMRDBottomWall((*h))) {
    fNSMRDBottomUnused++;
  }
  else if (COMET::IGeomInfo::SMRD().IsSMRDLeftWall((*h))) {
    fNSMRDLeftUnused++;
  }			
  else if (COMET::IGeomInfo::SMRD().IsSMRDRightWall((*h))) {
    fNSMRDRightUnused++;
  }	
      }
    }
  }


  
}


//*****************************************************
void COMET::IGlobalReconModule::FillFirstLastHits(COMET::IHitSelection& hits, IGlobalPID& globalObject){
//*****************************************************


  // sort hits in increasing Z
  std::sort(hits.begin(), hits.end(), SortHitsInZ);
    
  IGlobalHit* hitSaved[10];
  bool firstX=false;
  bool lastX=false;
  bool firstY=false;
  bool lastY=false;



  for(COMET::IHitSelection::iterator hit = hits.begin(); hit != hits.end(); hit++) {
    if (globalObject.NHitsSaved>=10 || (firstX && firstY)) break;

    if ((*hit)->IsXHit() && !firstX){
      hitSaved[globalObject.NHitsSaved] = new((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) IGlobalHit;    
      FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]);     
      globalObject.NHitsSaved++;
      firstX=true;
      if((*hit)->IsYHit()) firstY=true;
    }
    else if ((*hit)->IsYHit() && !firstY){
      hitSaved[globalObject.NHitsSaved] = new((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) IGlobalHit;    
      FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]); 
      globalObject.NHitsSaved++;
      firstY=true;
      if((*hit)->IsXHit()) firstX=true;
    }
  }

  for(COMET::IHitSelection::reverse_iterator hit = hits.rbegin(); hit != hits.rend(); hit++) {
    if (globalObject.NHitsSaved>=10 || (lastX && lastY)) break;

    if ((*hit)->IsXHit() && !lastX){
      hitSaved[globalObject.NHitsSaved] = new((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) IGlobalHit;    
      FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]);     
      globalObject.NHitsSaved++;
      lastX=true;
      if((*hit)->IsYHit()) lastY=true;
    }
    else if ((*hit)->IsYHit() && !lastY){
      hitSaved[globalObject.NHitsSaved] = new((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) IGlobalHit;    
      FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]); 
      globalObject.NHitsSaved++;
      lastY=true;
      if((*hit)->IsXHit()) lastX=true;
    }
  }

  

}

//*****************************************************
void COMET::IGlobalReconModule::FillOutermostHits(COMET::IHitSelection& hits, double charge_cut, IOutermostHits& outer){
//*****************************************************


  double minX =  1000000; 
  double maxX = -1000000; 
  double minY =  1000000; 
  double maxY = -1000000; 
  double minZ =  1000000; 
  double maxZ = -1000000; 

  COMET::IHandle<COMET::IHit> hitSave[6];
  

  for(COMET::IHitSelection::iterator hit = hits.begin(); hit != hits.end(); hit++) {
    
    if ((*hit)->GetCharge()<charge_cut) continue;
    
    if ((*hit)->IsXHit()){
      if ((*hit)->GetPosition().X() > maxX){
  maxX = (*hit)->GetPosition().X();
  hitSave[0] = *hit;
      }
      if ((*hit)->GetPosition().X() < minX){
  minX = (*hit)->GetPosition().X();
  hitSave[1] = *hit;
      }
    }
    if ((*hit)->IsYHit()){
      if ((*hit)->GetPosition().Y() > maxY){
  maxY = (*hit)->GetPosition().Y();
  hitSave[2] = *hit;
      }
      if ((*hit)->GetPosition().Y() < minY){
  minY = (*hit)->GetPosition().Y();
  hitSave[3] = *hit;
      }
    }
    if ((*hit)->IsZHit()){
      if ((*hit)->GetPosition().Z() > maxZ){
  maxZ = (*hit)->GetPosition().Z();
  hitSave[4] = *hit;
      }
      if ((*hit)->GetPosition().Z() < minZ){
  minZ = (*hit)->GetPosition().Z();
  hitSave[5] = *hit;
      }
    }    
  }

  FillGlobalHit(hitSave[0], outer.hitMaxX); 
  FillGlobalHit(hitSave[1], outer.hitMinX); 
  FillGlobalHit(hitSave[2], outer.hitMaxY); 
  FillGlobalHit(hitSave[3], outer.hitMinY); 
  FillGlobalHit(hitSave[4], outer.hitMaxZ); 
  FillGlobalHit(hitSave[5], outer.hitMinZ); 

  
}


//*****************************************************
void COMET::IGlobalReconModule::FillDelayedFgdClusters(COMET::ICOMETEvent& event, COMET::IHandle<COMET::IReconObjectContainer> allObjects){
//*****************************************************
  
  COMET::IHandle<COMET::IAlgorithmResult> ReconTracker = event.GetFit("ReconTracker");
  
  // __________________________________________________________
  // Fill the delayed FGD clusters.
  COMET::IHandle<COMET::IHitSelection> fgdHits;
  // Get all FGD hits
  if (ReconTracker)
    fgdHits = ReconTracker->GetHitSelection("FGD:unused");
  // Loop over vertices.
  if (fgdHits) {
    std::vector<COMET::IHitSelection> fgd_clusters;    

    // Find the earliest track time.
    double earliest_track_time = DEFAULT_MIN;

    if(allObjects->size() == 0)                                     // No tracks, so save all FGD clusters.
      earliest_track_time = 0.0;
    else {
      for(COMET::IReconObjectContainer::iterator track = allObjects->begin();
    track != allObjects->end(); track++) {
  double track_time = MedianObjectTime(*track);	

  if(track_time < earliest_track_time)
    earliest_track_time = track_time;
      }
    }    

    fEarliestTrackMedianHitTime = earliest_track_time;
    //if there are no rec track with FGD hits, set the track time to the mean vertex time
    if (earliest_track_time == 9999) fEarliestTrackMedianHitTime = 80;

    // Need to make this parameters somewhere else!!!
    
    double delay_time = 200*unit::ns;
    //double delay_time = 100*unit::ns;
    
    double min_distance = 5*unit::cm;
    double min_charge = 5.0;
    size_t min_hits = 2;
    for(COMET::IHitSelection::iterator hit = fgdHits->begin();
  hit != fgdHits->end(); hit++) {
      

      
      // Find all delayed hits.
      double dtime = (*hit)->GetTime() - earliest_track_time;
      
      // Need to make this parameter!!!
      if(dtime > delay_time) {	
  // Check if this hit is close to an existing cluster
  // of FGD hits.
  bool found_match = false;
  
  for(unsigned int i = 0; i < fgd_clusters.size(); i++) {
    for(COMET::IHitSelection::iterator chit = fgd_clusters[i].begin();
        chit != fgd_clusters[i].end(); chit++) {
      if(((*hit)->IsXHit() && (*chit)->IsXHit()) ||
         ((*hit)->IsYHit() && (*chit)->IsYHit())) {
        double pos_diff = ((*hit)->GetPosition() - (*chit)->GetPosition()).Mag();
        if(pos_diff < min_distance && (*hit)->GetCharge() > min_charge)
    found_match = true;
      }
    }
    
    if(found_match) {
      fgd_clusters[i].push_back(*hit);
      break;
    }
  }
  
  // If we didn't find a match, then start new hit selection.
  if(!found_match) {
    IHitSelection hits;
    hits.push_back(*hit);
    fgd_clusters.push_back(hits);
  }
      }
    }
    


    //Remove from the cluster the hits with a time value outside the range
    // (cluster average time - 2 * time rms, cluster average time + 2 * time rms)

    //loop over all the clusters
    for(unsigned int i = 0; i < fgd_clusters.size(); i++) {
      double avg_time_1 = 0.0;
      double sumT_1 = 0.0;
      double rmsTime_1 = 0.0; 

      //loop over the hits of the cluster
      for(COMET::IHitSelection::iterator hit = fgd_clusters[i].begin();
    hit != fgd_clusters[i].end(); hit++) { 
  //compute the average time of the cluster 
  avg_time_1 += (*hit)->GetTime();
      }
      avg_time_1 *= 1.0/((double)fgd_clusters[i].size());
      
      //loop over the hits of the cluster
      for(COMET::IHitSelection::iterator hit = fgd_clusters[i].begin();
    hit != fgd_clusters[i].end(); hit++) { 
  //compute the time rms of the cluster
  double squaredTime = pow((*hit)->GetTime() - avg_time_1, 2);
  sumT_1 += squaredTime;
      }
      sumT_1 *= 1.0/((double)fgd_clusters[i].size() - 1);
      rmsTime_1 = sqrt(sumT_1);
      
      //erase the "wrong" hits
      std::vector<COMET::IHitSelection::iterator> hits_to_erase;      
      for(COMET::IHitSelection::iterator hit = fgd_clusters[i].begin();
    hit != fgd_clusters[i].end(); hit++) {	
  double hitTime = (*hit)->GetTime();
  double minT = avg_time_1 - (2 * rmsTime_1);
  double maxT = avg_time_1 + (2 * rmsTime_1);
  if ((hitTime>maxT) || (hitTime<minT)) 
    hits_to_erase.push_back(hit);
      }    
      for (std::vector<COMET::IHitSelection::iterator>::iterator h = hits_to_erase.begin();
     h != hits_to_erase.end(); h++){
  fgd_clusters[i].erase(*h);
      } 
    }
    
  
   

    // For every delayed FGD cluster with more than one hit,
    // create a new TFgdClusterObject
    for(unsigned int i = 0; i < fgd_clusters.size(); i++) {
      if(fgd_clusters[i].size() > min_hits) {
  // Calculate average position and time.
  TVector3 avg_pos;
  double avg_time = 0.0;
  double total_charge = 0.0;
  TLorentzVector position_variance(1.0e8,1.0e8,0.0,0.0);
  
  double hitMinTime = 9999.0;           
  double hitMaxTime = -9999.0;          
  double hitMaxDistance = -9999.0;      
  std::vector<double> maxDistance;      

  for(COMET::IHitSelection::iterator hit = fgd_clusters[i].begin();
      hit != fgd_clusters[i].end(); hit++) {
    avg_pos += (*hit)->GetPosition();
    avg_time += (*hit)->GetTime();
    total_charge += (*hit)->GetCharge();
    if((*hit)->IsXHit())
      position_variance[0] = 0.0;
    else if((*hit)->IsYHit())
      position_variance[1] = 0.0;
    
    double hitTime = (*hit)->GetTime();                  
    if (hitTime < hitMinTime) hitMinTime = hitTime;      
    if (hitTime > hitMaxTime) hitMaxTime = hitTime;      

    
    TVector3 pos = (*hit)->GetPosition();
    double maxDist = -9999.;
    
    for(COMET::IHitSelection::iterator otherHit = fgd_clusters[i].begin();
        otherHit != fgd_clusters[i].end(); otherHit++) {
      TVector3 otherPos = (*(otherHit))->GetPosition();	    
      double dist = (pos - otherPos).Mag();
      if (dist > maxDist)
        maxDist = dist;
    }
    maxDistance.push_back(maxDist);	
  }

  for(unsigned int jj=0; jj<maxDistance.size(); jj++) {
    if (maxDistance[jj] > hitMaxDistance)
      hitMaxDistance = maxDistance[jj];
  }

  
  IFgdCluster* fgdCluster =
    new((*fDelayedClusters)[fNDelayedClusters++]) IFgdCluster;
  
  avg_pos *= 1.0/((double)fgd_clusters[i].size());
  avg_time *= 1.0/((double)fgd_clusters[i].size());
  
  double sumT = 0.0;                    
  double sumX = 0.0;                    
  double sumY = 0.0;          
  double sumZ = 0.0;          

  for(COMET::IHitSelection::iterator hit = fgd_clusters[i].begin();
      hit != fgd_clusters[i].end(); hit++) {
    double squaredTime = pow((*hit)->GetTime() - avg_time, 2);         
    sumT += squaredTime;                                               
    double squaredX = pow((*hit)->GetPosition().X() - avg_pos.X(), 2); 
    double squaredY = pow((*hit)->GetPosition().Y() - avg_pos.Y(), 2); 
    double squaredZ = pow((*hit)->GetPosition().Z() - avg_pos.Z(), 2); 
    sumX += squaredX;                                                  
    sumY += squaredY;                                                  
    sumZ += squaredZ;                                                  
  }
    
  sumT *= 1.0/((double)fgd_clusters[i].size() - 1);
  sumX *= 1.0/((double)fgd_clusters[i].size() - 1);
  sumY *= 1.0/((double)fgd_clusters[i].size() - 1);
  sumZ *= 1.0/((double)fgd_clusters[i].size() - 1);
  
  fgdCluster->TotalCharge = total_charge;
  fgdCluster->Position = TLorentzVector(avg_pos.X(),avg_pos.Y(),avg_pos.Z(),avg_time);
  fgdCluster->Variance = position_variance;
  
  fgdCluster->TimeDisp = hitMaxTime - hitMinTime; 
  fgdCluster->SpatialDisp = hitMaxDistance;       
  fgdCluster->NHits = fgd_clusters[i].size();     
  
  fgdCluster->PosRMS = TLorentzVector(sqrt(sumX), sqrt(sumY), sqrt(sumZ), sqrt(sumT)); 
  
      }
    }
  }
}



//*****************************************************
void COMET::IGlobalReconModule::FillFgdTimeBins(COMET::ICOMETEvent& event){
//*****************************************************

  bool isMC = event.GetContext().IsMC();

  fNFgdTimeBins = fgdUtils::getNumberTimeBins(event);
  for (int ibin=0; ibin<fNFgdTimeBins; ibin++) {

    IFgdTimeBin* bin = new((*fFgdTimeBins)[ibin]) IFgdTimeBin;
    COMET::IHandle<COMET::IAlgorithmResult> tbResult = fgdUtils::getTimeBinResult(event,ibin);

    bin->minTime = 2.e10;
    bin->maxTime = -2.e10;

    for (int ifgd=0; ifgd<2; ifgd++) {
      bin->nHits[ifgd] = 0;
      bin->rawChargeSum[ifgd] = 0.;
      bin->chargeWeightedPos[ifgd] = TVector3(0.,0.,0.);
      for (int ilayer=0; ilayer<30; ilayer++) {
  bin->chargePerLayer[ifgd][ilayer] = 0.;
      }
    }

    bin->NFGD1Unused=0;
    bin->NFGD2Unused=0;

    bin->FGD1Unused->Clear();
    bin->FGD2Unused->Clear();

    if (tbResult) {
      COMET::IHandle<COMET::IRealDatum> binStartTime = tbResult->Get<COMET::IRealDatum>("binStartTime");
      if (binStartTime) {
  bin->minTime = binStartTime->GetValue();
      }

      COMET::IHandle<COMET::IRealDatum> binEndTime = tbResult->Get<COMET::IRealDatum>("binEndTime");
      if (binEndTime) {
  bin->maxTime = binEndTime->GetValue();
      }

      // Fill the FGD unused hits
      COMET::IHandle<COMET::IHitSelection> fgdUnused = tbResult->GetHitSelection("unused");
      if (fgdUnused){	
  for(COMET::IHitSelection::iterator hit = fgdUnused->begin(); hit != fgdUnused->end(); hit++) {
    IGlobalHit* gHit = NULL;
    if      (COMET::IGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()))
      gHit = new((*bin->FGD1Unused)[bin->NFGD1Unused++]) IGlobalHit;
    else if (COMET::IGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))
      gHit = new((*bin->FGD2Unused)[bin->NFGD2Unused++]) IGlobalHit;
    else continue;

    FillGlobalHit(*hit, *gHit);
    
  }
      }

      // Also include hits from 2D fgdIsoRecon tracks for now
      COMET::IHandle<COMET::IReconObjectContainer> xz = tbResult->GetResultsContainer("xzUnmatchedFgdReconTracks");
      if (xz) {
  for(COMET::IReconObjectContainer::iterator tt = xz->begin();
      tt != xz->end(); tt++) {
    COMET::IHandle<COMET::IReconBase> object = *tt;
    if (object) {
      COMET::IHandle<COMET::IHitSelection> hits = object->GetHits();
      if (hits) {
        for(COMET::IHitSelection::iterator hit = hits->begin(); hit != hits->end(); hit++) {
    IGlobalHit* gHit = NULL;
    if      (COMET::IGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()))
      gHit = new((*bin->FGD1Unused)[bin->NFGD1Unused++]) IGlobalHit;
    else if (COMET::IGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))
      gHit = new((*bin->FGD2Unused)[bin->NFGD2Unused++]) IGlobalHit;
    else continue;

    FillGlobalHit(*hit, *gHit);

        }
      }
    }
  }
      }
      COMET::IHandle<COMET::IReconObjectContainer> yz = tbResult->GetResultsContainer("yzUnmatchedFgdReconTracks");
      if (yz) {
  for(COMET::IReconObjectContainer::iterator tt = yz->begin();
      tt != yz->end(); tt++) {
    COMET::IHandle<COMET::IReconBase> object = *tt;
    if (object) {
      COMET::IHandle<COMET::IHitSelection> hits = object->GetHits();
      if (hits) {
        for(COMET::IHitSelection::iterator hit = hits->begin(); hit != hits->end(); hit++) {
    IGlobalHit* gHit = NULL;
    if      (COMET::IGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()))
      gHit = new((*bin->FGD1Unused)[bin->NFGD1Unused++]) IGlobalHit;
    else if (COMET::IGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))
      gHit = new((*bin->FGD2Unused)[bin->NFGD2Unused++]) IGlobalHit;
    else continue;

    FillGlobalHit(*hit, *gHit);

        }
      }
    }
  }
      }

      // Cumulative result for all time bins
      fNFGD1Unused += bin->NFGD1Unused;
      fNFGD2Unused += bin->NFGD2Unused;


      COMET::IHandle<COMET::IHitSelection> tbHits = tbResult->GetHitSelection("fgd");
      if (tbHits) {

  // Extract FGD1 and FGD2 hits
  COMET::IHitSelection hits[2];

  TVector3 chargeSum[2];
  for(COMET::IHitSelection::iterator hit = tbHits->begin(); hit != tbHits->end(); hit++) {
    double charge = (*hit)->GetCharge();
    TVector3 pos = (*hit)->GetPosition();
    TVector3 isHit((*hit)->IsXHit(),(*hit)->IsYHit(),(*hit)->IsZHit());
    //	  std::cout << "isHit = (" << isHit.X() << "," << isHit.Y() << "," << isHit.Z() << ")" << std::endl;
    bool isInFgd[2];
    isInFgd[0] = COMET::IGeomInfo::FGD().IsInFGD1((*hit)->GetPosition());
    isInFgd[1] = COMET::IGeomInfo::FGD().IsInFGD2((*hit)->GetPosition());
    for (int ifgd=0; ifgd<2; ifgd++) {
      //	    std::cout << "isInFgd[" << ifgd << "] = " << isInFgd[ifgd] << std::endl;
      if (isInFgd[ifgd]) {
        hits[ifgd].push_back(*hit);
        bin->rawChargeSum[ifgd] += (*hit)->GetCharge();

        // Fill chargePerLayer
        COMET::IGeometryId geomId = (*hit)->GetGeomId();
        int xory = COMET::IGeomInfo::FGD().GetXorY(geomId);
        int module = COMET::IGeomInfo::FGD().GetModule(geomId);
        int layer = module*2 + xory;
        bin->chargePerLayer[ifgd][layer] += (*hit)->GetCharge();
        //std::cout << "Z, XorY, Module, Layer = " << pos.Z() << ", " << xory << ", " << module << ", " << layer << std::endl;

        for (int ix=0; ix<3; ix++) {
    if (isHit(ix)) {
      //		  std::cout << "filling hit for component " << ix << " in FGD" << ifgd+1 << std::endl;
      bin->chargeWeightedPos[ifgd](ix) += charge*pos(ix);
      chargeSum[ifgd](ix) += charge;
    }
        }
      }
    }

  }

  for (int ifgd=0; ifgd<2; ifgd++) {
//	  for (int ilayer=0; ilayer<30; ilayer++) {
//	    std::cout << "In FGD " << ifgd+1 << ", time bin " << ibin << ", charge per layer[" << ilayer << "] = " << bin->chargePerLayer[ifgd][ilayer] << std::endl;
//	  }
    for (int ix=0; ix<3; ix++) {
      if (chargeSum[ifgd](ix) > 1.e-3) {
        bin->chargeWeightedPos[ifgd](ix) *= 1./chargeSum[ifgd](ix);
      } else {
        bin->chargeWeightedPos[ifgd](ix) = -1.e-10;
        //	      std::cout << "No charge in ifgd,ix = (" << ifgd << "," << ix << "), so not filling chargeWeightedPos" << std::endl;
      }
    }
    //	  std::cout << "chargeWeightedPos[" << ifgd << "] = (" << bin->chargeWeightedPos[ifgd].X() << "," << bin->chargeWeightedPos[ifgd].Y() << "," << bin->chargeWeightedPos[ifgd].Z() << ")" << std::endl;
    bin->nHits[ifgd] = hits[ifgd].size();
  }


  // Fill outermost hits for this time bin
  double charge_cut = 0;
  FillOutermostHits(hits[0], charge_cut, bin->FGD1OutermostHits);
  FillOutermostHits(hits[1], charge_cut, bin->FGD2OutermostHits);

        if(isMC)
          bin->g4ID = TrackTruthInfo::GetG4TrajIDHits(*tbHits);
        else
          bin->g4ID = 0;
      }
    }

  }

}



// Helper method
// Get the median FGD hit time for this track.
//*************************************************************
double COMET::IGlobalReconModule::MedianObjectTime(COMET::IHandle<COMET::IReconBase> object){
//*************************************************************
  std::vector<double> times;
  COMET::IHandle<COMET::IHitSelection> hits = object->GetHits();
  if(hits){
    // make a vector of all the times for the 'good' hits
    for (COMET::IHitSelection::iterator hit = (hits)->begin(); hit != (hits)->end(); hit++) {
      if ((COMET::IGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()) ||
     COMET::IGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))) {
      
  times.push_back((*hit)->GetTime());
      }      
    }
  }  
  std::sort(times.begin(), times.end());
  
  unsigned len = times.size();
  
  if (len==0)                                                           // There must be an
    //return 0.0;
    return DEFAULT_MIN;
  else if (len % 2 != 0) {
    return times[(len-1) / 2];
  }
  else {
    return (times[len/2] + times[len/2 - 1]) / 2.0;
  }
  COMETDebug(" end median track ");
}

//*************************************************************
void COMET::IGlobalReconModule::GetFGDSimpleVA( COMET::ICOMETEvent& event, COMET::IHandle<COMET::IReconBase>& object, TLorentzVector& vertexPos){
//*************************************************************

  //time bin loop - match reco. time bin to track to get associated fgd hit selection
  int fgdRecoBinIndex = -1;
  if( !fgdUtils::GetTrackRecoTimeBinIndex( event, object, fgdRecoBinIndex ) )
    return;

  COMET::IHandle<COMET::IAlgorithmResult> binRecon = fgdUtils::getTimeBinResult(event,fgdRecoBinIndex);
  if( !binRecon )
    return;

  COMET::IHandle<COMET::IHitSelection> fgdHits = binRecon->GetHitSelection("fgd");
  if( !fgdHits )
    return;

  //run vertex activity measurement on track
  fgdUtils::SimpleVertexActivity(0, fgdHits, object, vertexPos );

  return;
}

//*************************************************************
void COMET::IGlobalReconModule::FillTPCPID(COMET::IHandle<COMET::IReconBase> object) {
//*************************************************************


  int dets = GetDetectorNumber(object);

  // create a new IGlobalPID
  ITpcPID* tpcObject = new((*fTPCPIDs)[fNTPCPIDs++]) ITpcPID;


  // Fill the true particle when it is found
  double pur, eff;
  COMET::IHandle<COMET::IG4Trajectory> G4track= GetG4Trajectory(*object,pur,eff);
  FillTrueParticle(G4track, pur, eff, tpcObject->TrueParticle);

  // extrapolate to the vertex (True vertex for the moment !!!!)
  tpcObject->PositionAtTrueVertex  = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
  tpcObject->PositionVarAtTrueVertex  = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
  tpcObject->DirectionAtTrueVertex  = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
  /*
  if (G4track && (object->UsesDetector(COMET::IReconBase::kTPC)  || object->UsesDetector(COMET::IReconBase::kFGD))){
    COMET::IHandle<COMET::IReconState> tstate = object->GetState();  
    if (tstate){
      bool ok = true;
      TVector3 snorm(0,0,1);
      ok = COMET::tman().PropagateToSurface(*tstate,G4track->GetInitialPosition().Vect(),snorm,*tstate);
      if (ok){
  tpcObject->Position  = TrackingUtils::GetPosition(tstate);
  tpcObject->Variance  = TrackingUtils::GetPositionVariance(tstate);     
  tpcObject->Direction = TrackingUtils::GetDirection(tstate);     
    }
    }
  }
  */
  //------- copy information from IReconPID to IGlobalPID --------
  
  tpcObject->AlgorithmName = object->GetAlgorithmName();
  tpcObject->Status     = object->CheckStatus(object->kSuccess);
  tpcObject->Chi2       = object->GetQuality();
  tpcObject->NDOF       = object->GetNDOF();
  tpcObject->NNodes     = object->GetNodes().size();
  tpcObject->Detectors  = dets;

  tpcObject->NHits = DEFAULT_MAX;  
  if (object->GetHits())
    tpcObject->NHits    = object->GetHits()->size();
  

  // The number of constituents. For composite objects
  tpcObject->NConstituents = 0;
  if (object->GetConstituents())
    tpcObject->NConstituents  = object->GetConstituents()->size();
  

  tpcObject->isForward     = (TrackingUtils::GetDirection(object).Z() > 0);  
  tpcObject->MomentumAtTrueVertex      = TrackingUtils::GetMomentum(object);
  tpcObject->MomentumErrorAtTrueVertex = TrackingUtils::GetMomentumError(object);
  tpcObject->Charge        = TrackingUtils::GetCharge(object);	
  tpcObject->EDeposit      = TrackingUtils::GetEDeposit(object);
  tpcObject->Width         = TrackingUtils::GetWidth(object).X(); //! take only the first component
  tpcObject->Cone          = TrackingUtils::GetCone(object);
  tpcObject->Length        = DEFAULT_MAX;
    
  
  // Set the PID back and front position, variance, direction and momentum
  COMET::IHandle<COMET::IReconState> firstState = TrackingUtils::GetFirstState(*object);
  if (firstState){
    tpcObject->FrontPosition       = TrackingUtils::GetPosition(firstState);
    tpcObject->FrontPositionVar    = TrackingUtils::GetPositionVariance(firstState);     
    tpcObject->FrontDirection      = TrackingUtils::GetDirection(firstState);     
    tpcObject->FrontDirectionVar   = TrackingUtils::GetDirectionVariance(firstState);     
    tpcObject->FrontMomentum       = TrackingUtils::GetMomentum(firstState);
    tpcObject->FrontMomentumError  = TrackingUtils::GetMomentumError(firstState);
  }
  else{
    tpcObject->FrontPosition       = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->FrontPositionVar    = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->FrontDirection      = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->FrontDirectionVar   = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->FrontMomentum       = DEFAULT_MAX; 
    tpcObject->FrontMomentumError  = DEFAULT_MAX; 
  }
  
  COMET::IHandle<COMET::IReconState> lastState  = TrackingUtils::GetLastState(*object);
  if (lastState){
    tpcObject->BackPosition       = TrackingUtils::GetPosition(lastState);
    tpcObject->BackPositionVar    = TrackingUtils::GetPositionVariance(lastState);     
    tpcObject->BackDirection      = TrackingUtils::GetDirection(lastState);     
    tpcObject->BackDirectionVar   = TrackingUtils::GetDirectionVariance(lastState);     
    tpcObject->BackMomentum       = TrackingUtils::GetMomentum(lastState);
    tpcObject->BackMomentumError  = TrackingUtils::GetMomentumError(lastState);
  }
  else{
    tpcObject->BackPosition       = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->BackPositionVar    = TLorentzVector(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->BackDirection      = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->BackDirectionVar   = TVector3(DEFAULT_MAX,DEFAULT_MAX,DEFAULT_MAX);
    tpcObject->BackMomentum       = DEFAULT_MAX; 
    tpcObject->BackMomentumError  = DEFAULT_MAX; 
  }
  
  // Set the PID information
  COMET::IHandle<COMET::IReconPID> PID = object;
  if (PID){      
    tpcObject->Likelihoods.push_back(PID->GetPIDWeight());
    COMET::IReconObjectContainer::const_iterator it;
    for (it=PID->GetAlternates().begin(); it!=PID->GetAlternates().end(); it++){
      COMET::IHandle<COMET::IReconPID> alter = *it;
      tpcObject->Likelihoods.push_back(alter->GetPIDWeight());
    }
  }
    
}


//*****************************************************
COMET::IHandle<COMET::IG4Trajectory> COMET::IGlobalReconModule::GetG4Trajectory(const COMET::IReconBase &object, double& purity, double& eff){
//*****************************************************

  purity = -0xABCDEF;
  eff = -0xABCDEF;

  //  std::cout<<"****** GET ASSOCIATE TRAJECTORY TO THE TRACK ****** "<<std::endl;
  COMET::IG4Trajectory traj;
  COMET::IHandle<COMET::IHitSelection> hits = object.GetHits();
  if (!hits ){
    return COMET::IHandle<COMET::IG4Trajectory>();
  }
  if(hits->size()<1){
    return COMET::IHandle<COMET::IG4Trajectory>();
  }

  int G4TrkId = TrackTruthInfo::GetG4TrajIDHits(*hits, eff, purity);
  const COMET::ICOMETEvent& event = *(COMET::IEventFolder::GetCurrentEvent());
  COMET::IHandle<COMET::IG4TrajectoryContainer> G4trajectories = event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");
  if (!G4trajectories){
    return COMET::IHandle<COMET::IG4Trajectory>();
  }
  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 ){
      COMET::IHandle<COMET::IG4Trajectory> traj(new COMET::IG4Trajectory((*trajectory)));
        return traj;
      }
  }
  

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

//*****************************************************
void COMET::IGlobalReconModule::UpdateCoincidences(COMET::IHandle<COMET::IMCHit> mch, COMET::IMCDigit* mcd, std::vector<int>& coinc, int& nhits){
//*****************************************************

  std::vector<IG4VHit*>  contrib;
  if (mch) 
    contrib = mch->GetContributors();
  else if (mcd)
    contrib = mcd->GetContributors();
  else
    return;
  
  
//	std::cout << "loop over hit contribuitors" << std::endl;
  for (std::vector<COMET::IG4VHit*>::const_iterator h = contrib.begin(); h != contrib.end(); ++h) {
    if( !*h ) continue; 
    COMET::IG4HitSegment *g4hitseg = dynamic_cast<COMET::IG4HitSegment*> (*h);
    if( g4hitseg != NULL ) {
      for( unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
  unsigned int vv = g4hitseg->GetContributors()[i];
  if( vv >= coinc.size() ) 
    coinc.resize(vv+1);
  coinc[vv]++;
  nhits++;
      }
    }    
  }
}

//*****************************************************
COMET::IHandle<COMET::IG4Trajectory> COMET::IGlobalReconModule::GetMainContributor(COMET::IHandle<COMET::IMCHit> mch, COMET::IHandle<COMET::IG4TrajectoryContainer> trajectories){
//*****************************************************

  if (!mch) return COMET::IHandle<COMET::IG4Trajectory>();
  if (!trajectories) return COMET::IHandle<COMET::IG4Trajectory>();

  std::vector<int> coinc; 
  coinc.resize(trajectories->size());

//	std::cout << "loop over hit contribuitors" << std::endl;
  for (std::vector<COMET::IG4VHit*>::const_iterator h = (mch->GetContributors()).begin(); h != (mch->GetContributors()).end(); ++h) {
    if( !*h ) continue; 
    COMET::IG4HitSegment *g4hitseg = dynamic_cast<COMET::IG4HitSegment*> (*h);
    if( g4hitseg != NULL ) {
      for( unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
  unsigned int vv = g4hitseg->GetContributors()[i];
  if( vv >= coinc.size() ) 
    coinc.resize(vv+1);
  coinc[vv]++;
      }
    }    
  }


  int maxim = 0; 
  int imax = -1;
  for( unsigned int i = 0; i < coinc.size() ; i++ ){
    if( coinc[i] > maxim ) {
      maxim = coinc[i];
      imax = i; 
    }
  }
  
  if( imax >= 0 ){
    COMET::IHandle<COMET::IG4Trajectory> traj(new COMET::IG4Trajectory((*trajectories)[imax])); 
    return traj;
  }

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

}


//*****************************************************
void COMET::IGlobalReconModule::GetHitsAssociatedToG4Trajectory(const COMET::IHitSelection& hits, COMET::IHandle<COMET::IG4Trajectory> traj, COMET::IHitSelection& traj_hits){
//*****************************************************
  
  if (!traj) return;

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

  if (!trajectories) return;

  // Loop over reconstructed hits 
  for( COMET::IHitSelection::const_iterator ht = hits.begin(); ht != hits.end(); ht++ ) { 
    
    COMET::IHandle<COMET::IComboHit> ch(*ht);
    if(ch){
      // Loop over hits in Combo
      for (COMET::IHitSelection::const_iterator ch_member = ch->GetHits().begin();      
     ch_member != ch->GetHits().end(); ++ch_member){
  
  // Get monte Carlo hit associated. 
  COMET::IHandle<COMET::IMCHit> mch;	  	 
  COMET::IHandle<COMET::IMultiHit> mhit = *ch_member;
  if (mhit)
    mch = *(mhit->begin());
  else
    mch = *ch_member;	      
  

  COMET::IHandle<COMET::IG4Trajectory> mainCont = GetMainContributor(mch,trajectories);
  if (mainCont)
    if (GetPointer(mainCont) == GetPointer(traj))
      traj_hits.push_back(mch);
      }
    }
    else{
      COMET::IHandle<COMET::IMultiHit> mhit = *ht;
      COMET::IHandle<COMET::IMCHit> mch;
      
      if (mhit)
  mch = *(mhit->begin());
      else
  mch = *ht;

  COMET::IHandle<COMET::IG4Trajectory> mainCont = GetMainContributor(mch,trajectories);
  if (mainCont)
    if (GetPointer(mainCont) == GetPointer(traj))
      traj_hits.push_back(mch);

    }
  }
  
}  

//*****************************************************
std::string COMET::IGlobalReconModule::GetObjectType(COMET::IHandle<COMET::IReconBase> object){
//*****************************************************

  std::string reco_class_name = "";

  COMET::IHandle<COMET::IReconObjectContainer> conContainer = object->GetConstituents();
  if (conContainer){
    COMET::IHandle<COMET::IReconBase> object2 = *(conContainer->begin());
    if (object2)             
      reco_class_name = object2->ClassName();
    
  }

  return reco_class_name;

}


//*****************************************************
bool COMET::IGlobalReconModule::IsTrackLike(COMET::IHandle<COMET::IReconBase> object){
//*****************************************************

  COMET::IHandle<COMET::IReconObjectContainer> conContainer = object->GetConstituents();  
  if (!conContainer) return true;
  if (conContainer->size()>1) return true;

  COMET::IHandle<COMET::IReconBase> object2 = *(conContainer->begin());
  if (!object2) return true;             

  std::string reco_class_name = object2->ClassName();
  
  if (reco_class_name.find("COMET::IReconPID")   !=std::string::npos || 
      reco_class_name.find("COMET::IReconTrack") !=std::string::npos )
    return true;
  else
    return false;

}


//*****************************************************************************
void  COMET::IGlobalReconModule::GetConstituentsInTracker(COMET::IHandle<COMET::IReconBase> t1, COMET::IReconObjectContainer& trackerObjects){
//*****************************************************************************
  
  // check if this object is tracker only. If so add it 
  if(IsTrackerOnly(t1)){
    trackerObjects.push_back(t1);
    return;
  }

  COMET::IHandle<COMET::IReconObjectContainer> t1const = t1->GetConstituents();
  if(!t1const) return;

  // loop over constituents
  COMET::IReconObjectContainer::const_iterator it1;
  for(it1 = t1const->begin(); it1 != t1const->end(); it1 ++){
    COMET::IHandle<COMET::IReconBase> object = *it1;
    if(!object) continue;
    GetConstituentsInTracker(object,trackerObjects);
  }							 
}

//*****************************************************************************
bool COMET::IGlobalReconModule::IsTrackerOnly(COMET::IHandle<COMET::IReconBase> t1){
//*****************************************************************************
  
  if(t1->UsesDetector(COMET::IReconBase::kP0D) || 
     t1->UsesDetector(COMET::IReconBase::kDSECal) || 
     t1->UsesDetector(COMET::IReconBase::kECal) || 
     t1->UsesDetector(COMET::IReconBase::kSMRD)){
    return false;
  }

  return true;
}

//*****************************************************************************
void COMET::IGlobalReconModule::DoAssociationBetweenTrackerAndGlobalObjects(const COMET::IReconObjectContainer& globalObjects, const COMET::IReconObjectContainer& trackerObjects){
//*****************************************************************************


  // 1. loop over all global objects
  for(COMET::IReconObjectContainer::const_iterator tt = globalObjects.begin();
      tt != globalObjects.end(); tt++){

    COMET::IHandle<COMET::IReconBase> gobject = *tt;
    if(!gobject) continue;
    
    // 2. Get the tracker objects in the global object
    COMET::IReconObjectContainer trackerObjects2;
    GetConstituentsInTracker(gobject,trackerObjects2);
    
    // 3. loop over tracker objects in the global object
    for(COMET::IReconObjectContainer::iterator tt2 = trackerObjects2.begin();
  tt2 != trackerObjects2.end(); tt2++){ 

      COMET::IHandle<COMET::IReconBase> tobject2 = *tt2;
      if(!tobject2) continue;
      
      // 4. loop over objects from ReconTracker
      for(COMET::IReconObjectContainer::const_iterator tt1 = trackerObjects.begin(); tt1 != trackerObjects.end(); tt1++){ 

  COMET::IHandle<COMET::IReconBase> tobject1 = *tt1;
  if(!tobject1) continue;

  // 5. Associate the tracker object in the global object with the ReconTracker object
  if(fabs(tobject1->GetQuality() - tobject2->GetQuality()) < 1e-6 && tobject1->GetDetectors() == tobject2->GetDetectors()){
    // this is the link between a ReconTracker object and a global object index
    fTrackerGlobalIndexMap[tobject1] = fGlobalIndexMap[gobject];
    break;
  }
      }
    }
  }
}

//*****************************************************************************
COMET::IHandle<COMET::IReconBase> COMET::IGlobalReconModule::GetTrackerReconVersionOfFGDIsoTrack(COMET::IHandle<COMET::IReconBase> object) {
//*****************************************************************************

  COMET::IHandle<COMET::IReconBase> trackerIsoObject;
  if (!object->UsesDetector(COMET::IReconBase::kFGD)||
      object->UsesDetector(COMET::IReconBase::kTPC)||
      object->UsesDetector(COMET::IReconBase::kECal)||
      object->UsesDetector(COMET::IReconBase::kDSECal)||
      object->UsesDetector(COMET::IReconBase::kP0D)||
      object->UsesDetector(COMET::IReconBase::kSMRD)
      ){
    return trackerIsoObject;
  }

  bool printStuff = false;
  double objecttime = -1.e10;
  COMET::IHandle<COMET::IReconObjectContainer> cons = object->GetConstituents();
  int ncons = cons->size();
  if (printStuff) {
    std::cout << "Found an FGD-only track" << std::endl;
    if (cons) {
      std::cout << "  and has " << ncons << " constituents" << std::endl;
    } else {
      std::cout << "  and has no constituents" << std::endl;
    }
  }
  COMET::IHandle<COMET::IReconTrack> track = object;
  COMET::IHandle<COMET::IReconPID> pid = object;
  if (track) {
    if (printStuff) std::cout << "   object is a track" << std::endl;
    objecttime = track->GetPosition().T();
  } else if (pid) {
    if (printStuff) std::cout << "   object is a PID" << std::endl;
    objecttime = pid->GetPosition().T();
  } else {
    std::cout << "Not a track or a PID?" << std::endl;
  }
  if (printStuff) {
    std::cout << "object time is " << objecttime << std::endl;
    for (int i=0; i<ncons; i++) {
      std::cout << "cons " << i << " detectors: P"
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kP0D)
    << " T"
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kTPC1)
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kFGD1)
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kTPC2)
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kFGD2)
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kTPC3)
    << " E"
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kECal)
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kDSECal)
    << " S"
    << (*cons)[i]->UsesDetector(COMET::IReconBase::kSMRD)
    << std::endl;
    }
  }

  // Find corresponding PID in ReconTracker "final" container
  COMET::IHandle<COMET::IAlgorithmResult> ReconTracker = COMET::IEventFolder::GetCurrentEvent()->GetFit("ReconTracker");
  if (ReconTracker){

    // Get the results from reconstruction
    COMET::IHandle<COMET::IReconObjectContainer> trackerObjects
      = ReconTracker->GetResultsContainer("final");

    if (trackerObjects) {
      for(COMET::IReconObjectContainer::iterator tt = trackerObjects->begin();
    tt != trackerObjects->end(); tt++) {
  if ((*tt)->UsesDetector(COMET::IReconBase::kFGD)&&
      !(*tt)->UsesDetector(COMET::IReconBase::kTPC)){
    COMET::IHandle<COMET::IReconTrack> track = *tt;
    COMET::IHandle<COMET::IReconPID> pid = *tt;
    double fgdtime = -1.e10;
    if (track) {
      std::cout << "   fgd object is a track?!?!?!?" << std::endl;
      fgdtime = track->GetPosition().T();
    } else if (pid) {
      if (printStuff) std::cout << "   fgd object is a PID" << std::endl;
      fgdtime = pid->GetPosition().T();
    } else {
      std::cout << "fgd not a track or a PID?" << std::endl;
    }
    if (printStuff) std::cout << "fgd object time is " << fgdtime << std::endl;
    if (fabs(fgdtime-objecttime) < 1.e-3){
      if (printStuff) std::cout << "replacing..." << std::endl;
      trackerIsoObject = *tt;
      break;
    }
  } else {
    if (printStuff) std::cout << "---no TPC component" << std::endl;
  }
      }
    } else {
      std::cout << "Why doesn't ReconTracker final results container exist?!?!?!?!?" << std::endl;
    }
  } else {
    std::cout << "Why doesn't ReconTracker result exist?!?!?!?!? " << std::endl;
  }
  if (!trackerIsoObject && printStuff) {
    std::cout << "!!!!!!!!!!!!!!!! ERROR:  FGD-only object NOT found!!!!!!!!!!!!!!!" << std::endl;
  }

  return trackerIsoObject;

}


//*****************************************************************************
void COMET::IGlobalReconModule::FillBrokenTracksMap(const COMET::IReconObjectContainer& globalObjects){
//*****************************************************************************

  //This function find the right links between broken tracks

  // clear the map from the previous event
  fBrokenIndexMap.clear();

  // map of uniqueIDs for the actual broken tracks
  std::map < COMET::IHandle<COMET::IReconBase>, std::vector<UInt_t> > brokenMap1;

  // map of uniqueIDs for the tracks associated to a broken track
  std::map < COMET::IHandle<COMET::IReconBase>, std::vector<UInt_t> > brokenMap2;

  // loop over global objects
  for(COMET::IReconObjectContainer::const_iterator tt = globalObjects.begin();
      tt != globalObjects.end(); tt++){

    COMET::IHandle<COMET::IReconBase> object = *tt;
    if(!object) continue;

    // get the broken IDs recursively
    std::vector<UInt_t> brokenIDs1;
    std::vector<UInt_t> brokenIDs2;
    bool found = GetBrokenIDs(object,brokenIDs1,brokenIDs2);

    // if this track is broken fill the maps
    if (found){
      brokenMap1[object]=brokenIDs1;
      brokenMap2[object]=brokenIDs2;
      
      if (debug_broken){
  std:: cout << " - " << object->GetUniqueID() << " " << object << std::endl;
  std:: cout << "     --> broken IDs1: ";
  for (unsigned int i=0;i<brokenIDs1.size();i++){
    std::cout << " " << brokenIDs1[i];
  }
  std:: cout << std::endl;
  std:: cout << "     --> broken IDs2: ";
  for (unsigned int i=0;i<brokenIDs2.size();i++){
    std::cout << " " << brokenIDs2[i];
  }
  std:: cout << std::endl;
      }
    }

  }

  // now do the actual association
  std::map< COMET::IHandle<COMET::IReconBase>, std::vector<UInt_t> >::const_iterator tt1;
  std::map< COMET::IHandle<COMET::IReconBase>, std::vector<UInt_t> >::const_iterator tt2;
  for(tt1 = brokenMap1.begin();  tt1 != brokenMap1.end(); tt1++){
    
    COMET::IHandle<COMET::IReconBase> object1 = tt1->first;
    if(!object1) continue;

    const std::vector<UInt_t>& brokenIDs1 = tt1->second;

    std::vector<UInt_t> brokenIDsFinal;
    bool found=false;
    // does any of the brokenIDs of this object exist as a independent global object ?
    for(tt2 = brokenMap2.begin();  tt2 != brokenMap2.end(); tt2++){
      
      COMET::IHandle<COMET::IReconBase> object2 = tt2->first;
      if(!object2) continue;

      const std::vector<UInt_t>& brokenIDs2 = tt2->second;
      
      for (unsigned int i=0;i<brokenIDs1.size();i++){
  for (unsigned int j=0;j<brokenIDs2.size();j++){
    if (brokenIDs1[i] == brokenIDs2[j]){
      brokenIDsFinal.push_back(object2->GetUniqueID());
      found=true;
      break;
    }
  }
  if (found) break;
      }
      if (found) break;
    }
    
    // fill the map with top level association 
    if (found){
      fBrokenIndexMap[object1]=brokenIDsFinal;    

      if (debug_broken){      
  std:: cout << " - " << object1->GetUniqueID() << " " << object1 << std::endl;
  std:: cout << "     --> broken IDs FINAL: ";
  for (unsigned int i=0;i<brokenIDsFinal.size();i++){
    std::cout << " " << brokenIDsFinal[i];
  }
  std:: cout << std::endl;
      }
    }
  }
}

//*****************************************************************************
bool COMET::IGlobalReconModule::GetBrokenIDs(COMET::IHandle<COMET::IReconBase> object, std::vector<UInt_t>& brokenIDs1, std::vector<UInt_t>& brokenIDs2){
//*****************************************************************************

  bool found=false;
  
  // get the broken IDs for this object
  COMET::IHandle<COMET::IDataVector> broken =  object->Get< COMET::IDataVector >("brokenID");
  if (broken){
    brokenIDs1.push_back(object->GetUniqueID());
    found=true;
    for (unsigned int i=0;i<broken->size();i++){
      UInt_t ID =  broken->At<COMET::IIntegerDatum>(i)->GetValue(); 
      brokenIDs2.push_back(ID);
    }
  }
  
  // Get The BrokenIDs for the original object
  if(ReconObjectUtils::GetOriginalObject(object)){
    found = GetBrokenIDs(ReconObjectUtils::GetOriginalObject(object),brokenIDs1,brokenIDs2);
  }
  else{
    // get the broken IDs for the constituents (only if the track has no original object)
    COMET::IHandle<COMET::IReconObjectContainer> cons = object->GetConstituents();
    if (!cons) return found;
    
    for(COMET::IReconObjectContainer::const_iterator tt = cons->begin();
  tt != cons->end(); tt++){
      
      // only for PID constituents
      COMET::IHandle<COMET::IReconPID> object2 = *tt;
      if(!object2) continue;
      
      bool found2 = GetBrokenIDs(object2,brokenIDs1,brokenIDs2);
      if (found2) 
  found=true;
    }
  }

  if (found) 
    brokenIDs1.push_back(object->GetUniqueID());
  
  return found;

}