#include "IDatum.hxx"
#include "IRealDatum.hxx"
#include "TPrincipal.h"
#include "IGeomInfo.hxx"
#include "TrackTruthInfo.hxx"
#include "IShowerTruthInfo.hxx"

#include <sstream>

#include "ITrackerReconECalModule.hxx"

ClassImp(COMET::ITrackerReconECalModule);
ClassImp(COMET::ITrackerReconECalModule::IReconECalTrack);
ClassImp(COMET::ITrackerReconECalModule::IReconECalShower);
ClassImp(COMET::ITrackerReconECalModule::IReconECalObject);
ClassImp(COMET::ITrackerReconECalModule::IReconECalCluster);
ClassImp(COMET::ITrackerReconECalModule::IReconECalUnmatchedObject);

#define CVSTAG "\
  $Name: v5r19 $"
#define CVSID  "\
  $Id: ITrackerReconECalModule.cxx,v 1.79 2012/11/05 11:48:25 marks Exp $"

COMET::ITrackerReconECalModule::ITrackerReconECalModule(const char *name,
const char *title) {

  SetNameTitle(name, title);
  // Enable this module by default:
  fIsEnabled = kTRUE;
  fDescription = "Tracker ECAL recon output";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;
  fIsMC = false;
  total = 0;
  // Tree variables
  fReconObjects = new TClonesArray("COMET::ITrackerReconECalModule::IReconECalObject", 200);
  fUnmatchedObjects = new TClonesArray("COMET::ITrackerReconECalModule::IReconECalUnmatchedObject", 200);
}


COMET::ITrackerReconECalModule::~ITrackerReconECalModule() 
{

}

Bool_t COMET::ITrackerReconECalModule::ProcessFirstEvent(COMET::ICOMETEvent& event) {
  return true;
}

bool COMET::ITrackerReconECalModule::IsFullEventWorthSaving(COMET::ICOMETEvent& event){
  return (fNReconObjects ? true : false);
}

void COMET::ITrackerReconECalModule::InitializeModule() {

}

void COMET::ITrackerReconECalModule::InitializeBranches() {

  fOutputTree->Branch("NReconObject", &fNReconObjects, "NReconObjects/I", fBufferSize);
  fOutputTree->Branch("NReconTrackLike", &fNReconTrackLike, "NReconTrackLike/I", fBufferSize);
  fOutputTree->Branch("NReconShowerLike", &fNReconShowerLike, "NReconShowerLike/I", fBufferSize);
  fOutputTree->Branch("NReconCluster", &fNReconCluster, "NReconCluster/I", fBufferSize);
  fOutputTree->Branch("NUnmatchedObject", &fNUnmatchedObjects, "NUnmatchedObjects/I", fBufferSize);
  // A TClonesArray of IReconECalObject objects
  fOutputTree->Branch("ReconObject", &fReconObjects, fBufferSize, fSplitLevel);
  fOutputTree->Branch("UnmatchedObject", &fUnmatchedObjects, fBufferSize, fSplitLevel);

}

bool COMET::ITrackerReconECalModule::FillTree(COMET::ICOMETEvent& event) {

  fIsMC = event.GetContext().IsMC();

  // Initialise
  fNReconObjects = 0;
  fNUnmatchedObjects = 0;
  fNReconTrackLike = 0;
  fNReconShowerLike = 0;
  fNReconCluster = 0;

  fReconObjects->Clear();
  fUnmatchedObjects->Clear();

  // Getting the final Algorithm Result from ReconECal.
  COMET::IHandle<COMET::IAlgorithmResult> ecalFinalResult = event.GetFit("ReconECal");


  if(ecalFinalResult) {
  for(COMET::IDataVector::iterator dataIter = ecalFinalResult->begin();
  dataIter != ecalFinalResult->end(); ++dataIter) {
    std::string class_name = (*dataIter)->ClassName();
    TString object_name = (*dataIter)->GetName();

    // Looking for ObjectContainers that are labelled correctly
    if(class_name.find("COMET::IReconObjectContainer") != std::string::npos && object_name.Contains("final")) {
      COMET::IReconObjectContainer* reconContainer =
        dynamic_cast<COMET::IReconObjectContainer*>(*dataIter);

      if(!reconContainer)
        COMETWarn("Dynamic cast of ECALPid ReconObjectContainer didn't work");

      // Getting bunch number.
      std::string std_object_name = (*dataIter)->GetName();
      std::string::size_type idx = std_object_name.find_last_of('_');
      std::string bunch_id = std_object_name.substr(idx+1);

      // Loop looking for IReconPID's in the top recon conatiner
      for(COMET::IReconObjectContainer::iterator reconIter =
        reconContainer->begin();
      reconIter != reconContainer->end();  ++reconIter) {

        std::string reco_class_name = (*reconIter)->ClassName();
  
  // Get IReconPID's
  if(reco_class_name.find("COMET::IReconPID") != std::string::npos) {
    
          COMET::IHandle<COMET::IReconPID> pid = (*reconIter);

    COMET::IHandle<COMET::IReconObjectContainer> conContainer = pid->GetConstituents();
    for(COMET::IReconObjectContainer::iterator conIter =
    conContainer->begin();
        conIter != conContainer->end();  ++conIter) { 

      std::string constituent_class_name = (*conIter)->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 = (*conIter);
        IReconECalObject *ReconECalObject = 0;
        
        if(track->GetEDeposit()) {
    // Making a IECALObject ClonesArray entry...
    ReconECalObject = new((*fReconObjects)[fNReconObjects]) IReconECalObject;
    //
    // Get the cluster information
    //
    COMET::IHandle<COMET::IReconObjectContainer> conClusterContainer = track->GetConstituents();
    for(COMET::IReconObjectContainer::iterator conClusterIter = conClusterContainer->begin();
        conClusterIter != conClusterContainer->end();  ++conClusterIter) { 

      std::string cluster_class_name = (*conClusterIter)->ClassName();

      COMET::IHandle<COMET::IReconCluster> cluster = (*conClusterIter);
      ReconECalObject->Cluster.Status = cluster->GetStatus();
      ReconECalObject->Cluster.NDOF = cluster->GetNDOF();
      ReconECalObject->Cluster.Quality = cluster->GetQuality();
      ReconECalObject->Cluster.Position = cluster->GetPosition();
      ReconECalObject->Cluster.PositionVar = cluster->GetPositionVariance();
      ReconECalObject->Cluster.EDeposit = cluster->GetEDeposit();

      fNReconCluster++;
    }
    //
    // These bools are to show what the ReconECal using default 
    // cuts thinks the object is, a shower, or a track.
    ReconECalObject->IsShowerLike = false;
    ReconECalObject->IsTrackLike = true;

    ReconECalObject->TimeBunch = atoi(bunch_id.c_str());
    
    // If tracks have hits grab some hit lvl information.  It may
    // make sence in the future to remove some of this info when
    // we have a working hit lvl information module.

                // Save the unique object ID to allow matching to global recon objects
                ReconECalObject->UniqueID = pid->GetUniqueID();

    // Calculate the pointing vector.
    TVector3 Pointing(-9999, -9999, -9999);
    if(track->GetHits()){
      COMET::IHandle<COMET::IHitSelection> hitsel = track->GetHits();

                  // Get the 2D and 3D radon results.
                  ReconECalObject->FillRadonInfo(pid);

      ReconECalObject->FillHitInfo(hitsel);
      Pointing = ReconECalObject->CalculatePointing(hitsel);
                  if(fIsMC){
                      ReconECalObject->G4ID = TrackTruthInfo::GetG4TrajIDHits(*hitsel);

                      // Fill the truth-matching information from the new truth-matching algorithms.
                      COMET::IShowerTruthInfo showerTruth(pid);

                      std::vector<COMET::IShowerTruthInfo::TruthInfo> results = showerTruth.GetPrimaryTrajectories();
                      if (results.size() > 0) {
                        ReconECalObject->G4ID_Primary = results[0].ParticleID;
                        ReconECalObject->Completeness_Primary = results[0].Completeness_NumHits;
                        ReconECalObject->Cleanliness_Primary = results[0].Cleanliness_NumHits;
                      }

                      results = showerTruth.GetRecursiveTrajectories();
                      if (results.size() > 0) {
                        ReconECalObject->G4ID_Recursive = results[0].ParticleID;
                        ReconECalObject->Completeness_Recursive = results[0].Completeness_NumHits;
                        ReconECalObject->Cleanliness_Recursive = results[0].Cleanliness_NumHits;
                      }
                      if (results.size() > 1) {
                        ReconECalObject->G4ID_Recursive2 = results[1].ParticleID;
                      }
                      if (results.size() > 2) {
                        ReconECalObject->G4ID_Recursive3 = results[2].ParticleID;
                      }
                      if (results.size() > 3) {
                        ReconECalObject->G4ID_Recursive4 = results[3].ParticleID;
                      }

                      results = showerTruth.GetSingleTrajectories();
                      if (results.size() > 0) {
                        ReconECalObject->G4ID_Single = results[0].ParticleID;
                        ReconECalObject->Completeness_Single = results[0].Completeness_NumHits;
                        ReconECalObject->Cleanliness_Single = results[0].Cleanliness_NumHits;
                      }
                  }
    }
    ReconECalObject->Pointing = Pointing;

    // Collecting information stored in a IReconTrack..
    ReconECalObject->Track.Curvature = track->GetCurvature();
    ReconECalObject->Track.Direction = track->GetDirection();
    ReconECalObject->Track.EDeposit = track->GetEDeposit();
    ReconECalObject->Track.NDOF = track->GetNDOF();
    ReconECalObject->Track.Position = track->GetPosition();
    ReconECalObject->Track.PositionVar = track->GetPositionVariance();
    ReconECalObject->Track.Quality = track->GetQuality();
    ReconECalObject->Track.Width = track->GetWidth();
    ReconECalObject->Track.BackPosition = ReconECalObject->FindBackPosition(track->GetHits());
    ReconECalObject->Track.Status = track->GetStatus();

    fNReconObjects++;
                total++;
    fNReconTrackLike++;

    // Adding the Information stored in Datums in the Track
    // This for example includes parametres used in the Pid and 
    // Energy Recon. 
    ReconECalObject->FillDatumInfo(track);

    // Adding properties of the Alternative EMShower Hypothesis 
    for(COMET::IReconObjectContainer::const_iterator altIter =
          pid->GetAlternates().begin();
        altIter != pid->GetAlternates().end();  ++altIter) { 		  
      std::string alternates_class_name = (*altIter)->ClassName();

      if (alternates_class_name.find("COMET::IReconPID") != std::string::npos) {

        COMET::IHandle<COMET::IReconPID> altPid = (*altIter);
        COMET::IHandle<COMET::IReconObjectContainer> altConContainer = altPid->GetConstituents();	    
        for(COMET::IReconObjectContainer::iterator altConIter =
        altConContainer->begin();
      altConIter != altConContainer->end();  ++altConIter) { 
          
          std::string altConstituent_class_name = (*altConIter)->ClassName();
          
          if(altConstituent_class_name.find("COMET::IReconShower") != std::string::npos) {
      COMET::IHandle<COMET::IReconShower> shower = (*altConIter);
      if(shower->GetEDeposit()) {
        ReconECalObject->Shower.ConeAngle = shower->GetConeAngle();
        ReconECalObject->Shower.Direction = shower->GetDirection();
        ReconECalObject->Shower.EDeposit = shower->GetEDeposit();
        ReconECalObject->Shower.NDOF = shower->GetNDOF();
        ReconECalObject->Shower.Position = shower->GetPosition();
        ReconECalObject->Shower.PositionVar = shower->GetPositionVariance();
        ReconECalObject->Shower.Quality = shower->GetQuality();
        ReconECalObject->Shower.BackPosition = ReconECalObject->FindBackPosition(shower->GetHits());
        ReconECalObject->Shower.Status = shower->GetStatus();
    
      }//end if AltShowerEnergy =/= 0
          }// end if AltShower
        }// end for AltConContainer
      }// end if AltPid
    }// end for AltContainer
        }
      }
      // Now looking for TReconPid's that have Shower Contituents
      else if (constituent_class_name.find("COMET::IReconShower") != std::string::npos) {
        COMET::IHandle<COMET::IReconShower> shower = (*conIter);
        IReconECalObject *ReconECalObject = 0;

        if(shower->GetEDeposit()) {
    ReconECalObject = new((*fReconObjects)[fNReconObjects]) IReconECalObject;
    //
    // Get the cluster information
    //
    COMET::IHandle<COMET::IReconObjectContainer> conClusterContainer = shower->GetConstituents();
    for(COMET::IReconObjectContainer::iterator conClusterIter = conClusterContainer->begin();
        conClusterIter != conClusterContainer->end();  ++conClusterIter) { 

      std::string cluster_class_name = (*conClusterIter)->ClassName();

      COMET::IHandle<COMET::IReconCluster> cluster = (*conClusterIter);
      ReconECalObject->Cluster.Status = cluster->GetStatus();
      ReconECalObject->Cluster.NDOF = cluster->GetNDOF();
      ReconECalObject->Cluster.Quality = cluster->GetQuality();
      ReconECalObject->Cluster.Position = cluster->GetPosition();
      ReconECalObject->Cluster.PositionVar = cluster->GetPositionVariance();
      ReconECalObject->Cluster.EDeposit = cluster->GetEDeposit();

      fNReconCluster++;
    }
    //

    ReconECalObject->IsShowerLike = true;
    ReconECalObject->IsTrackLike = false;

    ReconECalObject->TimeBunch = atoi(bunch_id.c_str());

    // If showers has hits grab some hit lvl information.  It may
    // make sence in the future to remove some of this info when
    // we have a working hit lvl information module.

                // Save the unique object ID to allow matching to global recon objects
                ReconECalObject->UniqueID = pid->GetUniqueID();

    // Calculate the pointing vector.
                TVector3 Pointing(-9999, -9999, -9999);

    if(shower->GetHits()){
      COMET::IHandle<COMET::IHitSelection> hitsel = shower->GetHits();

                  // Get the 2D and 3D radon results.
                  ReconECalObject->FillRadonInfo(pid);

      ReconECalObject->FillHitInfo(hitsel);
                  Pointing = ReconECalObject->CalculatePointing(hitsel);
                  if(fIsMC){
                      ReconECalObject->G4ID = TrackTruthInfo::GetG4TrajIDHits(*hitsel);

                      // Fill the truth-matching information from the new truth-matching algorithms.
                      COMET::IShowerTruthInfo showerTruth(pid);

                      std::vector<COMET::IShowerTruthInfo::TruthInfo> results = showerTruth.GetPrimaryTrajectories();
                      if (results.size() > 0) {
                        ReconECalObject->G4ID_Primary = results[0].ParticleID;
                        ReconECalObject->Completeness_Primary = results[0].Completeness_NumHits;
                        ReconECalObject->Cleanliness_Primary = results[0].Cleanliness_NumHits;
                      }

                      results = showerTruth.GetRecursiveTrajectories();
                      if (results.size() > 0) {
                        ReconECalObject->G4ID_Recursive = results[0].ParticleID;
                        ReconECalObject->Completeness_Recursive = results[0].Completeness_NumHits;
                        ReconECalObject->Cleanliness_Recursive = results[0].Cleanliness_NumHits;
                      }
                      if (results.size() > 1) {
                        ReconECalObject->G4ID_Recursive2 = results[1].ParticleID;
                      }
                      if (results.size() > 2) {
                        ReconECalObject->G4ID_Recursive3 = results[2].ParticleID;
                      }
                      if (results.size() > 3) {
                        ReconECalObject->G4ID_Recursive4 = results[3].ParticleID;
                      }

                      results = showerTruth.GetSingleTrajectories();
                      if (results.size() > 0) {
                        ReconECalObject->G4ID_Single = results[0].ParticleID;
                        ReconECalObject->Completeness_Single = results[0].Completeness_NumHits;
                        ReconECalObject->Cleanliness_Single = results[0].Cleanliness_NumHits;
                      }


                  }
          }

                ReconECalObject->Pointing = Pointing;

    ReconECalObject->Shower.ConeAngle = shower->GetConeAngle();
    ReconECalObject->Shower.Direction = shower->GetDirection();
    ReconECalObject->Shower.EDeposit = shower->GetEDeposit();
    ReconECalObject->Shower.NDOF = shower->GetNDOF();
    ReconECalObject->Shower.Position = shower->GetPosition();
    ReconECalObject->Shower.PositionVar = shower->GetPositionVariance();
    ReconECalObject->Shower.Quality = shower->GetQuality();
    ReconECalObject->Shower.BackPosition = ReconECalObject->FindBackPosition(shower->GetHits());
    ReconECalObject->Shower.Status = shower->GetStatus();

    fNReconObjects++;
                total++;
    fNReconShowerLike++;
    
    // Adding the Information stored in Datums in the Shower
    // This for example includes parametres using in TMVA and Energy
    // fitting. 
    ReconECalObject->FillDatumInfo(shower);

    // Adding properties of the EMTrack Hypothesis 
    for(COMET::IReconObjectContainer::const_iterator altIter =
          pid->GetAlternates().begin();
        altIter != pid->GetAlternates().end();  ++altIter) { 		  
      std::string alternates_class_name = (*altIter)->ClassName();

      if (alternates_class_name.find("COMET::IReconPID") != std::string::npos) {

        COMET::IHandle<COMET::IReconPID> altPid = (*altIter);
        COMET::IHandle<COMET::IReconObjectContainer> altConContainer = altPid->GetConstituents();	    
        for(COMET::IReconObjectContainer::iterator altConIter =
        altConContainer->begin();
      altConIter != altConContainer->end();  ++altConIter) { 
          
          std::string altConstituent_class_name = (*altConIter)->ClassName();
          
          if(altConstituent_class_name.find("COMET::IReconTrack") != std::string::npos) {
      COMET::IHandle<COMET::IReconTrack> track = (*altConIter);
      if(track->GetEDeposit()) {
        
        ReconECalObject->Track.Curvature = track->GetCurvature();
        ReconECalObject->Track.Direction = track->GetDirection();
        ReconECalObject->Track.EDeposit = track->GetEDeposit();
        ReconECalObject->Track.NDOF = track->GetNDOF();
        ReconECalObject->Track.Position = track->GetPosition();
        ReconECalObject->Track.PositionVar = track->GetPositionVariance();
        ReconECalObject->Track.Quality = track->GetQuality();
        ReconECalObject->Track.Width = track->GetWidth();
        ReconECalObject->Track.BackPosition = ReconECalObject->FindBackPosition(track->GetHits());
        ReconECalObject->Track.Status = track->GetStatus();
       
      }//end if AltTrackEnergy =/= 0
          }// end if AltTrack
        }// end for AltConContainer
      }// end if AltPid
    }// end for AltContainer

        }
      }
    }    
  }
      }
    } // End section for "final"
    // Leigh: Now let's take care of the unmatched objects
    else if(class_name.find("COMET::IReconObjectContainer") != std::string::npos && object_name.Contains("other")) {
      COMET::IReconObjectContainer* unCont =	dynamic_cast<COMET::IReconObjectContainer*>(*dataIter);

      if(!unCont) COMETWarn("Dynamic cast of unmatched ReconObjectContainer didn't work");

      IReconECalUnmatchedObject *unmatchedObject = 0x0;

      for(COMET::IReconObjectContainer::iterator cl = unCont->begin(); cl != unCont->end(); ++cl){
        std::string unmatched_class_name = (*cl)->ClassName();
 
        // Look for COMET::IReconCluster objects 
        if(unmatched_class_name.find("COMET::IReconCluster") != std::string::npos) {
          COMET::IHandle<COMET::IReconCluster> clust = (*cl);
          if(clust->GetHits()){
            unmatchedObject = new((*fUnmatchedObjects)[fNUnmatchedObjects]) IReconECalUnmatchedObject;
            unmatchedObject->FillHitInfo(clust->GetHits());
            ++fNUnmatchedObjects;

                                    if(fIsMC){
                // Fill the truth-matching information from the new truth-matching algorithms.
                COMET::IShowerTruthInfo showerTruth(clust);
    
                std::vector<COMET::IShowerTruthInfo::TruthInfo> results = showerTruth.GetPrimaryTrajectories();
                if (results.size() > 0) {
                  unmatchedObject->G4ID_Primary = results[0].ParticleID;
                }
    
                results = showerTruth.GetRecursiveTrajectories();
                if (results.size() > 0) {
                  unmatchedObject->G4ID_Recursive = results[0].ParticleID;
                }
    
                results = showerTruth.GetSingleTrajectories();
                if (results.size() > 0) {
                  unmatchedObject->G4ID_Single = results[0].ParticleID;
                }
                                    }
          }
        }
      }
    } // End unmatched object section.
  }
  } // if(ecalFinalResult)
  return true;
}

COMET::ITrackerReconECalModule::IReconECalObject::IReconECalObject(){
  
    IsTrackLike = false;
    IsShowerLike = false;
    //IsEMLike = false;
    //IsHADLike = false;
    G4ID_Primary = -1;
    Completeness_Primary = -1;
    Cleanliness_Primary = -1;
    G4ID_Recursive = -1;
    Completeness_Recursive = -1;
    Cleanliness_Recursive = -1;
    G4ID_Single = -1;
    Completeness_Single = -1;
    Cleanliness_Single = -1;
}

void COMET::ITrackerReconECalModule::IReconECalObject::FillRadonInfo(COMET::IHandle<COMET::IReconPID> pid){

    // Get the radon line datums from the ecal PID object
    COMET::IHandle<COMET::IRealDatum> xy_first_radon_datum = pid->Get<COMET::IRealDatum>("XY_First_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> xy_second_radon_datum = pid->Get<COMET::IRealDatum>("XY_Second_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> xy_third_radon_datum = pid->Get<COMET::IRealDatum>("XY_Third_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> xz_first_radon_datum = pid->Get<COMET::IRealDatum>("XZ_First_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> xz_second_radon_datum = pid->Get<COMET::IRealDatum>("XZ_Second_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> xz_third_radon_datum = pid->Get<COMET::IRealDatum>("XZ_Third_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> yz_first_radon_datum = pid->Get<COMET::IRealDatum>("YZ_First_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> yz_second_radon_datum = pid->Get<COMET::IRealDatum>("YZ_Second_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> yz_third_radon_datum = pid->Get<COMET::IRealDatum>("YZ_Third_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> full_first_radon_datum = pid->Get<COMET::IRealDatum>("3D_First_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> full_second_radon_datum = pid->Get<COMET::IRealDatum>("3D_Second_Radon_Line");
    COMET::IHandle<COMET::IRealDatum> full_third_radon_datum = pid->Get<COMET::IRealDatum>("3D_Third_Radon_Line");

    // Copy the information into the correct oaAnalysis output variable
    this->Radon_XY_First_Prong_NHits = int(xy_first_radon_datum->GetVector()[0]); 
    this->Radon_XY_First_Prong_Angle = int(xy_first_radon_datum->GetVector()[1]); 
    this->Radon_XY_First_Prong_Distance = int(xy_first_radon_datum->GetVector()[2]); 
    this->Radon_XY_Second_Prong_NHits = int(xy_second_radon_datum->GetVector()[0]); 
    this->Radon_XY_Second_Prong_Angle = int(xy_second_radon_datum->GetVector()[1]); 
    this->Radon_XY_Second_Prong_Distance = int(xy_second_radon_datum->GetVector()[2]); 
    this->Radon_XY_Third_Prong_NHits = int(xy_third_radon_datum->GetVector()[0]); 
    this->Radon_XY_Third_Prong_Angle = int(xy_third_radon_datum->GetVector()[1]); 
    this->Radon_XY_Third_Prong_Distance = int(xy_third_radon_datum->GetVector()[2]); 

    this->Radon_XZ_First_Prong_NHits = int(xz_first_radon_datum->GetVector()[0]); 
    this->Radon_XZ_First_Prong_Angle = int(xz_first_radon_datum->GetVector()[1]); 
    this->Radon_XZ_First_Prong_Distance = int(xz_first_radon_datum->GetVector()[2]); 
    this->Radon_XZ_Second_Prong_NHits = int(xz_second_radon_datum->GetVector()[0]); 
    this->Radon_XZ_Second_Prong_Angle = int(xz_second_radon_datum->GetVector()[1]); 
    this->Radon_XZ_Second_Prong_Distance = int(xz_second_radon_datum->GetVector()[2]); 
    this->Radon_XZ_Third_Prong_NHits = int(xz_third_radon_datum->GetVector()[0]); 
    this->Radon_XZ_Third_Prong_Angle = int(xz_third_radon_datum->GetVector()[1]); 
    this->Radon_XZ_Third_Prong_Distance = int(xz_third_radon_datum->GetVector()[2]); 

    this->Radon_YZ_First_Prong_NHits = int(yz_first_radon_datum->GetVector()[0]); 
    this->Radon_YZ_First_Prong_Angle = int(yz_first_radon_datum->GetVector()[1]); 
    this->Radon_YZ_First_Prong_Distance = int(yz_first_radon_datum->GetVector()[2]); 
    this->Radon_YZ_Second_Prong_NHits = int(yz_second_radon_datum->GetVector()[0]); 
    this->Radon_YZ_Second_Prong_Angle = int(yz_second_radon_datum->GetVector()[1]); 
    this->Radon_YZ_Second_Prong_Distance = int(yz_second_radon_datum->GetVector()[2]); 
    this->Radon_YZ_Third_Prong_NHits = int(yz_third_radon_datum->GetVector()[0]); 
    this->Radon_YZ_Third_Prong_Angle = int(yz_third_radon_datum->GetVector()[1]); 
    this->Radon_YZ_Third_Prong_Distance = int(yz_third_radon_datum->GetVector()[2]); 

    TVector3 first_prong_direction(full_first_radon_datum->GetVector()[0], full_first_radon_datum->GetVector()[1], full_first_radon_datum->GetVector()[2]);
    this->Radon_3D_First_Prong_Direction = first_prong_direction;
    this->Radon_3D_First_Prong_NHits = int(full_first_radon_datum->GetVector()[3]);

    TVector3 second_prong_direction(full_second_radon_datum->GetVector()[0], full_second_radon_datum->GetVector()[1], full_second_radon_datum->GetVector()[2]);
    this->Radon_3D_Second_Prong_Direction = second_prong_direction;
    this->Radon_3D_Second_Prong_NHits = int(full_second_radon_datum->GetVector()[3]);

    TVector3 third_prong_direction(full_third_radon_datum->GetVector()[0], full_third_radon_datum->GetVector()[1], full_third_radon_datum->GetVector()[2]);
    this->Radon_3D_Third_Prong_Direction = third_prong_direction;
    this->Radon_3D_Third_Prong_NHits = int(full_third_radon_datum->GetVector()[3]);
}

void COMET::ITrackerReconECalModule::IReconECalObject::FillHitInfo(COMET::IHandle<COMET::IHitSelection> hitsel){

  // Setting hit lvl quantities to appropriate start values
  this->mostUpStreamLayerHit  = 36;
  this->mostDownStreamLayerHit   = -1;
  this->maxPerpBarHit = -1;
  this->maxParaBarHit = -1;
  this->minBarHit = 100;
  this->NLayersHit     = 0;
  this->NHits          = 0;
  this->TotalHitCharge = 0;
  this->AverageHitTime = 0;
  this->NIsXHits       = 0;
  this->NIsYHits       = 0;
  this->NIsZHits       = 0;
  this->ObjectLength = 0.0;
  this->AverageZPosition = 0.0;

  std::string module;

  /*
  // Added by Gavin Davies, BarHit and LayerHit as referred to in examineEcalHits.cxx
  this->BarHit         = 0;
  this->LayerHit       = 0;
  this->IsXHit         = 0;
  this->IsYHit         = 0;
  this->IsZHit         = 0;
  */

  bool layer[36] = {false};
  std::vector<TVector3> layerHitPos(36);
  std::vector<double> layerCharge(36,0.0);

  // Vector to store the layer for each hit
  std::vector<int> layerhits;
  //std::vector<int> layerhits(1000,-1);
  double HitCharge = -1;

  // Looping through hits and filling the Hit info
  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);
    Int_t barHit = COMET::IGeomInfo::Get().ECAL().GetBarNumber(reconHit->GetContributor(0)->GetGeomId());

    layerhits.push_back(layerHit);

    if(hit == hitsel->begin()){
        module = COMET::IGeomInfo::Get().ECAL().GetModule(reconHit->GetGeomId()).GetName();
    }

    // Average Z Position
    this->AverageZPosition += HitPos.Z();
    
    // Charge weighted layer positon
    TVector3 chargePos = reconHit->GetPosition() * reconHit->GetCharge();
    layerHitPos[layerHit] += chargePos;
    layerCharge[layerHit] += reconHit->GetCharge();

    /*
    // Added by Gavin, to fill IsXHit, IsYHit, IsZHit and LayerHit
    //this->IsXHit = reconHit->GetConstituents().at(0)->IsXHit();
    //this->IsYHit = reconHit->GetConstituents().at(0)->IsYHit();
    //this->IsZHit = reconHit->GetConstituents().at(0)->IsZHit();
    this->IsXHit = reconHit->IsXHit();
    this->IsYHit = reconHit->IsYHit();
    this->IsZHit = reconHit->IsZHit();
    this->LayerHit = layerHit;
    // Added to fill BarHit
    if(this->IsXHit){
    this->BarHit = COMET::IGeomInfo::ECAL().GetActiveXBar(HitPos);
        }
        if(this->IsYHit){
    this->BarHit = COMET::IGeomInfo::ECAL().GetActiveYBar(HitPos);;
        }

    */

    // Using the array to collect the number of layers hit.
    if( (!layer[layerHit]) ){
      layer[layerHit] = true;
      this->NLayersHit++;
    }
    // Filling numbers on hits in differenet types of bars.
    if(reconHit->GetContributor(0)->IsXHit()){
      this->NIsXHits++;
    }
    if(reconHit->GetContributor(0)->IsYHit()){
      this->NIsYHits++;
    }
    if(reconHit->GetContributor(0)->IsZHit()){
      this->NIsZHits++;
    }
    // Obtaining the max and min bars hit 
    if(barHit < this->minBarHit){
        this->minBarHit = barHit;
    }
    if(barHit > this->maxParaBarHit && !reconHit->GetContributor(0)->IsZHit()){
        this->maxParaBarHit = barHit;
    }
    if(barHit > this->maxPerpBarHit && reconHit->GetContributor(0)->IsZHit()){
        this->maxPerpBarHit = barHit;
    }
    // Obtaining the first and last layers hit.
    if(layerHit < this->mostUpStreamLayerHit){
      this->mostUpStreamLayerHit = layerHit;
    }
    if(layerHit > this->mostDownStreamLayerHit){
      this->mostDownStreamLayerHit = layerHit;
    }
    // Collecting other useful hit information
    this->NHits++;
    this->TotalHitCharge += reconHit->GetCharge();
    this->AverageHitTime += reconHit->GetTime() * reconHit->GetCharge();
    // Save the layer with the highest hit charge
    if(reconHit->GetCharge() > HitCharge){
      HitCharge = reconHit->GetCharge();
      this->MaxHitChargeLayer = layerHit+1;
    }
  }

  this->Module = module;

  this->AverageHitTime /= this->TotalHitCharge;
  this->AverageZPosition /= this->NHits;
  
  // Charge weighted layer position
  for(int i = 0; i < 36; ++i){
    if(layerCharge[i] != 0.0){
      layerHitPos[i] *= (1 / layerCharge[i]);
    }
  }
  
  // Now calculate object length
  TVector3 objLength = layerHitPos[this->mostDownStreamLayerHit] - layerHitPos[this->mostUpStreamLayerHit];
  this->ObjectLength = objLength.Mag();

  // Sort the layer for each hit by increasing order
  sort(layerhits.begin(),layerhits.end());
  // Temporary variables to be used
  int seclayer = -1;
  int doublelayer = -1;
  int doublelayersec = -1;
  int layertwohits = -1;
  int layerthreehits = -1;
  int layerfourhits = -1;
  int lastlayermanyhits = -1;
  int layerA = -1;
  int layerB = -1;
  int NLayers1Hits = -1;
  int NLayers2Hits = -1;
  int NLayers3Hits = -1;
  int NLayers4Hits = -1;
  int NLayers5Hits = -1;
  int NLayers6Hits = -1;
  int NLayers7Hits = -1;
  int maxhits = -1;

  for(unsigned int i=0; i < layerhits.size(); i++){
    // if(layerhits.at(i) == -1){continue;}
    int layer = layerhits.at(i);

    if(layer > layerA){
      layerA = layer;
      NLayers1Hits++;
      seclayer = -1;
    }
    else{
      seclayer++;
      doublelayer++;
      lastlayermanyhits = layer;
      if(layer > layerB){
        layerB = layer;
        doublelayersec++;
        if(doublelayersec == 0){
          layertwohits = layer;
        }
        else if(doublelayersec == 1){
          layerthreehits = layer;
        }
        else if(doublelayersec == 2){
          layerfourhits = layer;
        }
      }

      if(seclayer > maxhits)
  maxhits = seclayer;

      if(seclayer == 0)
        NLayers2Hits++;
      else if(seclayer == 1)
        NLayers3Hits++;
      else if(seclayer == 2)
        NLayers4Hits++;
      else if(seclayer == 3)
        NLayers5Hits++;
      else if(seclayer == 4)
        NLayers6Hits++;
      else if(seclayer == 5)
        NLayers7Hits++;
    }
  }

  // Save the variables - The first layer is always starting from 1
  this->FirstLayer = layerhits.front()+1;
  this->LastLayer = layerhits.back()+1;
  this->NHitsAtLayersWithManyHits = doublelayer+1;
  this->NLayersOneHits =  NLayers1Hits+1;
  this->NLayersTwoHits = NLayers2Hits+1;
  this->NLayersThreeHits = NLayers3Hits+1;
  this->NLayersFourHits = NLayers4Hits+1;
  this->NLayersFiveHits = NLayers5Hits+1;
  this->NLayersSixHits = NLayers6Hits+1;
  this->NLayersSevenHits = NLayers7Hits+1;
  this->FirstLayerManyHits = layertwohits+1;
  this->SecondLayerManyHits = layerthreehits+1;
  this->ThirdLayerManyHits = layerfourhits+1;
  this->LastLayerManyHits =  lastlayermanyhits+1;
  this->MaxHitsInALayer = maxhits+1;

  // Empty the vector
  layerhits.clear();

}

void COMET::ITrackerReconECalModule::IReconECalObject::FillDatumInfo(COMET::IHandle<COMET::IReconBase> base){

  // =========================Pid Vars==========================  
  //                         ----------
  // AMR..
  if(base->Get< COMET::IRealDatum >("AMR")) {
    this->PID_AMR = 
      base->Get< COMET::IRealDatum >("AMR")->GetValue();
  }
  else{
    this->PID_AMR = -1;
  }
  
  // Circularity
  // Combined and in each view
  if(base->Get< COMET::IRealDatum >("Circularity")) {
    this->PID_Circularity = 
      base->Get< COMET::IRealDatum >("Circularity")->GetValue();
  }
  else{
    this->PID_Circularity = -1;
  }
  if(base->Get< COMET::IRealDatum >("CircularityView0")) {
    this->PID_CircularityView0 = 
      base->Get< COMET::IRealDatum >("CircularityView0")->GetValue();
  }
  else{
    this->PID_CircularityView0 = -1;
  }
  if(base->Get< COMET::IRealDatum >("CircularityView1")) {
    this->PID_CircularityView1 = 
      base->Get< COMET::IRealDatum >("CircularityView1")->GetValue();
  }
  else{
    this->PID_CircularityView1 = -1;
  }
  
  // Containment
  // This isn't PID, but it's easy to save it the same way as a PID variable
  if(base->Get< COMET::IRealDatum >("Containment")) {
    this->Containment = 
      base->Get< COMET::IRealDatum >("Containment")->GetValue();
  }
  else{
    this->Containment = -1;
  }
  
  // Angle..
  if(base->Get< COMET::IRealDatum >("Angle")){ 
    this->PID_Angle = 
      base->Get< COMET::IRealDatum >("Angle")->GetValue();
  }
  else{
    this->PID_Angle = -1;
  }
  
  // Max Ratio..
  if(base->Get< COMET::IRealDatum >("Max_Ratio")){ 
    this->PID_Max_Ratio = 
      base->Get< COMET::IRealDatum >("Max_Ratio")->GetValue();
  }
  else {
    this->PID_Max_Ratio = -1;
  }
  
  // Truncated Max Ratio..
  if(base->Get< COMET::IRealDatum >("TruncatedMaxRatio")){
    this->PID_TruncatedMaxRatio =
      base->Get< COMET::IRealDatum >("TruncatedMaxRatio")->GetValue();
  }
  else {
    this->PID_TruncatedMaxRatio = -1;
  }

  // Transverse charge ratio.
  if(base->Get< COMET::IRealDatum >("TransverseChargeRatio")){
    this->PID_TransverseChargeRatio =
      base->Get< COMET::IRealDatum >("TransverseChargeRatio")->GetValue();
  }
  else {
    this->PID_TransverseChargeRatio = -1;
  }

  // NormalizedMipChi2.
  if(base->Get< COMET::IRealDatum >("NormalizedMipChi2")){
    this->PID_NormalizedMipChi2 =
      base->Get< COMET::IRealDatum >("NormalizedMipChi2")->GetValue();
  }
  else {
    this->PID_NormalizedMipChi2 = -1;
  }

  // NormalizedMipChi2AllLayers.
  if(base->Get< COMET::IRealDatum >("NormalizedMipChi2AllLayers")){
    this->PID_NormalizedMipChi2AllLayers =
      base->Get< COMET::IRealDatum >("NormalizedMipChi2AllLayers")->GetValue();
  }
  else {
    this->PID_NormalizedMipChi2AllLayers = -1;
  }

    //FrontBackRatio. Note - Currently assumes particle is coming from the tracker into the ECal
    if(base->Get< COMET::IRealDatum >("FrontBackRatio")){
    this->PID_FrontBackRatio =
      base->Get< COMET::IRealDatum >("FrontBackRatio")->GetValue();
  }
  else {
  COMET::IHandle<COMET::IHitSelection> hits = base->Get<COMET::IHitSelection>("hits");
  double FrontBackRatio(-1);
  double FrontBackRatioNumerator(-1);	
  double FrontBackRatioDenominator(-1);
  if(hits)
    {
      CalculateFrontBackRatio(hits, FrontBackRatio, FrontBackRatioNumerator, FrontBackRatioDenominator);
    }
  this->PID_FrontBackRatio = FrontBackRatio;
  this->PID_FrontBackRatioNumerator = FrontBackRatioNumerator;
  this->PID_FrontBackRatioDenominator = FrontBackRatioDenominator;
  }


  // Charge Ratio
  /* Seems not to be included currently, commenting out for now
  if(base->Get< COMET::IRealDatum >("ChargeRatio")){ 
    this->PID_ChargeRatio = 
      base->Get< COMET::IRealDatum >("ChargeRatio")->GetValue();
  }
  else{
    this->PID_ChargeRatio = -1;
  }
  */

  // Shower Angle
  if(base->Get< COMET::IRealDatum >("ShowerAngle")){ 
    this->PID_ShowerAngle = 
      base->Get< COMET::IRealDatum >("ShowerAngle")->GetValue();
  }
  else{
    this->PID_ShowerAngle = -1;
  }
  
  // Asymmetry
  if(base->Get< COMET::IRealDatum >("Asymmetry")){
    this->PID_Asymmetry = 
      base->Get< COMET::IRealDatum >("Asymmetry")->GetValue();
  }
  else{
    this->PID_Asymmetry = -1;
  }
  
  // Mean Position
  if(base->Get< COMET::IRealDatum >("MeanPosition")) {
    this->PID_MeanPosition = 
      base->Get< COMET::IRealDatum >("MeanPosition")->GetValue();
  }
  else{
    this->PID_MeanPosition = -1;
  }
  
  // Qskew
  if(base->Get< COMET::IRealDatum >("Qskew")){ 
    this->PID_Qskew = 
      base->Get< COMET::IRealDatum >("Qskew")->GetValue();
  }
  else {
    this->PID_Qskew = -1;
  }
  
  // Shower Width
  if(base->Get< COMET::IRealDatum >("ShowerWidth")){ 
    this->PID_ShowerWidth = 
      base->Get< COMET::IRealDatum >("ShowerWidth")->GetValue();
  }
  else{
    this->PID_ShowerWidth = -1;
  }
  
  // TMVA Output - value between 0 and 1 to choose a Track or Shower hypothesis
  if(base->Get< COMET::IRealDatum >("TrShval")){ 
    this->PID_TrShval = 
      base->Get< COMET::IRealDatum >("TrShval")->GetValue();
  }
  else{
    this->PID_TrShval = -1;
  }
  
  // Likelihood PID Output
  if(base->Get< COMET::IRealDatum >("LLR_MIP_EM")){
    this->PID_LLR_MIP_EM =
      base->Get< COMET::IRealDatum >("LLR_MIP_EM")->GetValue();
  }
  else{
    this->PID_LLR_MIP_EM = -9999.;
  }

  // Likelihood PID Output
  if(base->Get< COMET::IRealDatum >("LLR_MIP_Pion")){
    this->PID_LLR_MIP_Pion =
      base->Get< COMET::IRealDatum >("LLR_MIP_Pion")->GetValue();
  }
  else{
    this->PID_LLR_MIP_Pion = -9999.;
  }

  // Likelihood PID Output
  if(base->Get< COMET::IRealDatum >("LLR_EM_HIP")){
    this->PID_LLR_EM_HIP =
      base->Get< COMET::IRealDatum >("LLR_EM_HIP")->GetValue();
  }
  else{
    this->PID_LLR_EM_HIP = -9999.;
  }

  // Likelihood PID Output
  if(base->Get< COMET::IRealDatum >("LLR_MIP_EM_LowMomentum")){
    this->PID_LLR_MIP_EM_LowMomentum =
      base->Get< COMET::IRealDatum >("LLR_MIP_EM_LowMomentum")->GetValue();
  }
  else{
    this->PID_LLR_MIP_EM_LowMomentum = -9999.;
  }

  // EM Likelihood
  if(base->Get< COMET::IRealDatum >("EMLikelihood")){ 
    this->PID_EMLikelihood = 
      base->Get< COMET::IRealDatum >("EMLikelihood")->GetValue();
  }
  else{
    this->PID_EMLikelihood = -1;
  }
  
  // TMVA Output - value between 0 and 1 to choose a EM or Hadronic hypothesis
  if(base->Get< COMET::IRealDatum >("EMHadVal")) {
    this->PID_EMHadVal = base->Get< COMET::IRealDatum >("EMHadVal")->GetValue();
  }
  else{
    this->PID_EMHadVal = -1;
  }
  
  // Now Looking for TDatums in the Cluster Constituent of the Tracks and 
  // Showers.  These Datum include;
  //   >  Energy likelihood fit parametres, uncertainties etc...
  //   >  Michel Electron Tag information on any tagged delayed clusters  
  //   >  Thrust fit parameters
  //   >  Matching Likelihood from 3D Matching
  COMET::IHandle<COMET::IReconObjectContainer> baseContainer = base->GetConstituents();
  for(COMET::IReconObjectContainer::iterator baseIter =
  baseContainer->begin();
      baseIter != baseContainer->end();  ++baseIter) { 
    
    std::string constituent_class_name = (*baseIter)->ClassName();
    if(constituent_class_name.find("COMET::IReconCluster") != std::string::npos) {
      COMET::IHandle<COMET::IReconCluster> cluster = (*baseIter);
      
      // ====================Energy Variables========================
      //                    ------------------
  
      // EM Energy Fit Result   
      if(cluster->Get< COMET::IRealDatum >("EMEnergyFitResult")){
  this->EMEnergyFit_Result = cluster->Get< COMET::IRealDatum >("EMEnergyFitResult")->GetVector().front();
      }
      else {
  this->EMEnergyFit_Result = -1;
      }
      
      // EM Energy Fit Uncertainty
      if(cluster->Get< COMET::IRealDatum >("EMEnergyFitUncertainty")){
  this->EMEnergyFit_Uncertainty = cluster->Get< COMET::IRealDatum >("EMEnergyFitUncertainty")->GetValue();
      }
      else {
  this->EMEnergyFit_Uncertainty = -1;
      }
      
      // EM Energy Fit Likelihoods
      if(cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")){
  this->EMEnergyFit_Likelihood_energyGrad = cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")->GetVector()[0];
  this->EMEnergyFit_Likelihood_qsum_like = cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")->GetVector()[1];
  this->EMEnergyFit_Likelihood_energy_qsumGrad = cluster->Get< COMET::IRealDatum >("EMEnergyFitLikelihood")->GetVector()[2];
  
      }
      else{
  this->EMEnergyFit_Likelihood_energyGrad = -1;
  this->EMEnergyFit_Likelihood_energy_qsumGrad = -1;
  this->EMEnergyFit_Likelihood_qsum_like = -1;
      }

      // EM Energy Fit Parameters
      if(cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")){
  this->EMEnergyFit_Para_QSum = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[0];
  this->EMEnergyFit_Para_QMean = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[1];
  this->EMEnergyFit_Para_QRMS = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[2];
  this->EMEnergyFit_Para_QSkew = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[3];
  this->EMEnergyFit_Para_QMax = cluster->Get< COMET::IRealDatum >("EMEnergyFitParameters")->GetVector()[4];
      }
      else{
  this->EMEnergyFit_Para_QSum = -1;
  this->EMEnergyFit_Para_QMean = -1;
  this->EMEnergyFit_Para_QRMS = -1;
  this->EMEnergyFit_Para_QSkew = -1;
  this->EMEnergyFit_Para_QMax = -1;
      }


      // Michel Electron Tag Output
      if(cluster->Get< COMET::IRealDatum >("NDelayedClusters")){
  this->MElectronTag_NDelayedCluster = cluster->Get< COMET::IRealDatum >("NDelayedClusters")->GetVector()[0];
      }
      else{
  this->MElectronTag_NDelayedCluster= -1;
      }   
      if(MElectronTag_NDelayedCluster > 0){
  for(Int_t dclust = 0; dclust < 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 Variables from Datum

      MElectronTag_NDelayedHits.push_back(Delta[0]);
      MElectronTag_DeltaX.push_back(Delta[1]);
      MElectronTag_DeltaY.push_back(Delta[2]);
      MElectronTag_DeltaZ.push_back(Delta[3]);
      MElectronTag_DeltaT.push_back(Delta[4]);
      MElectronTag_EDep.push_back(Delta[5]);
      MElectronTag_TrackEnd.push_back(Delta[6]);
    }
  }
      }

      // Pid output from the Kalman Fitter
      if(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")){
  this->KFParameter = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[0];
  this->KFMultiTracksTPC  = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[1]);
  this->KFMultiTracksECAL = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[2]);
  this->KFNNodes = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[3]);
  this->KFParameterNodes = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[4];
  this->KFNDOF = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[5]);
  this->KFQuality = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[6];
  this->KFIsMatched = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[7];
  this->KFHitQuality = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatum")->GetVector()[8];
      }
      else{
  this->KFParameter     = -1;
  this->KFMultiTracksTPC  = -1;
  this->KFMultiTracksECAL = -1;
  this->KFNNodes = -1;
  this->KFParameterNodes = -1;
  this->KFQuality = -1;
  this->KFNDOF = -1;
  this->KFIsMatched = false;
  this->KFHitQuality = -1;
      }

      // Pid output from the Kalman Fitter (in ECAL Stand Alone mode
      if(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")){
  this->KFParameterEcalStandAlone = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[0];
  this->KFMultiTracksTPCEcalStandAlone  = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[1]);
  this->KFMultiTracksECALEcalStandAlone = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[2]);
  this->KFNNodesEcalStandAlone = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[3]);
  this->KFParameterNodesEcalStandAlone = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[4];
  this->KFNDOFEcalStandAlone = int(cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[5]);
  this->KFQualityEcalStandAlone = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[6];
  this->KFIsMatchedEcalStandAlone = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[7];
  this->KFHitQualityEcalStandAlone = cluster->Get< COMET::IRealDatum >("ECALTPCMatchingDatumEcalStandAlone")->GetVector()[8];
      }
      else{
  this->KFParameterEcalStandAlone     = -1;
  this->KFMultiTracksTPCEcalStandAlone  = -1;
  this->KFMultiTracksECALEcalStandAlone = -1;
  this->KFNNodesEcalStandAlone = -1;
  this->KFParameterNodesEcalStandAlone = -1;
  this->KFQualityEcalStandAlone = -1;
  this->KFNDOFEcalStandAlone = -1;
  this->KFIsMatchedEcalStandAlone = false;
  this->KFHitQualityEcalStandAlone = -1;
      }
     
      // ====================Thrust Variables======================== 
      //                    ------------------

      if(cluster->Get< COMET::IRealDatum >("Thrust")){
        std::vector<double> thrust = cluster->Get< COMET::IRealDatum >("Thrust")->GetVector();
  this->Thrust = thrust[0];
  this->ThrustOrigin.SetXYZ(thrust[1], thrust[2], thrust[3]);
  this->ThrustAxis.SetXYZ(thrust[4], thrust[5], thrust[6]);
      }
      else {
        this->Thrust = -1.;
      }

      // ====================Matching Likelihood======================== 
      //                    ------------------

      if(cluster->Get< COMET::IRealDatum >("Likelihood")){
        this->MatchingLikelihood = cluster->Get< COMET::IRealDatum >("Likelihood")->GetValue();
      }
      else {
        this->MatchingLikelihood = -1.;
      }



      break;  
    }
  }
}

TLorentzVector COMET::ITrackerReconECalModule::IReconECalObject::FindBackPosition(COMET::IHandle<COMET::IHitSelection> hits){

  TLorentzVector pos;
  double totalCharge(0);

  // Iterate over hits to get the total charge in the back layer.
  for(COMET::IHitSelection::iterator it = hits->begin(); it != hits->end(); ++it){
    COMET::IHandle<COMET::IReconHit> hit = *it;
    int layer = COMET::IGeomInfo::Get().ECAL().GetLayerNumber(hit->GetGeomId());

    if(layer == mostDownStreamLayerHit){
      totalCharge += hit->GetCharge();
    }
  }

  // Iterate over hits again to find a charge weighted position.
  for(COMET::IHitSelection::iterator it = hits->begin(); it != hits->end(); ++it){
    COMET::IHandle<COMET::IReconHit> hit = *it;
    int layer = COMET::IGeomInfo::Get().ECAL().GetLayerNumber(hit->GetGeomId());

    if(layer == mostDownStreamLayerHit){
    
      TLorentzVector temp(hit->GetPosition(),hit->GetTime());
      temp *= hit->GetCharge() / totalCharge;

      pos += temp;

    }
  }

//	std::cout << totalCharge << ", " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << std::endl;

  return pos;
}

TVector3 COMET::ITrackerReconECalModule::IReconECalObject::CalculatePointing(COMET::IHandle<COMET::IHitSelection> hits){
                                                         
  TPrincipal pca(3,"ND");

  // Add the charge weighted hits to the pca.
  for(COMET::IHitSelection::iterator hit = hits->begin(); hit != hits->end(); ++hit){
    double row[3] = { (*hit)->GetPosition().X(),
                      (*hit)->GetPosition().Y(),
                      (*hit)->GetPosition().Z() };
    for (double q=0; q < (*hit)->GetCharge(); ++q) {
      pca.AddRow(row);
    }
  }

  // Find the principal components of the shower.
  pca.MakePrincipals();

  // Find the centre of the shower.
  //   The P2X function does this:
  //   Get the mean position of the shower.
  //    for (Int_t i = 0; i < fNumberOfVariables; i++){
  //      x[i] = fMeanValues(i);
  //      for (Int_t j = 0; j < nTest; j++)
  //        x[i] += p[j] * (fIsNormalised ? fSigmas(i) : 1)
  //      * fEigenVectors(i,j);
  //    }
  //   Give p[j] = 0, so actually just return x[i], which is the
  //   mean value over all data points for x, y and z.
  double posP[] = {0.,0.,0.};
  double posX[3];
  pca.P2X(posP, posX, 3);
  TVector3 pcaCentre; pcaCentre.SetXYZ(posX[0], posX[1], posX[2]);
 
  // Set the primary pca axis.
  TVector3 pcaAxis; pcaAxis.SetXYZ((*pca.GetEigenVectors())[0][0], 
           (*pca.GetEigenVectors())[1][0], 
           (*pca.GetEigenVectors())[2][0]);
 
  // Use the PCA info to calcualte pointing.
  double sumcharge = 0;
  double numerator = 0;
  for(COMET::IHitSelection::iterator hit = hits->begin(); hit != hits->end(); ++hit){
    sumcharge += (*hit)->GetCharge();
    numerator += (*hit)->GetCharge() * pcaAxis.Dot((*hit)->GetPosition() - pcaCentre);
  }

  TVector3 pointing;
  pointing = (numerator / sumcharge) * pcaAxis;

  return pointing;
}



COMET::ITrackerReconECalModule::IReconECalUnmatchedObject::IReconECalUnmatchedObject() {
  G4ID_Primary = -1;
  G4ID_Recursive = -1;
  G4ID_Single = -1;
};

void COMET::ITrackerReconECalModule::IReconECalUnmatchedObject::FillHitInfo(COMET::IHandle<COMET::IHitSelection> hits)
{
  // Initialise the variables.
  this->NHits = 0;
  this->View = -1;
  this->TotalHitCharge = 0.0;
  this->AverageHitTime = 0.0;
  this->mostUpStreamLayerHit = 0;
  this->mostDownStreamLayerHit = 0;

  this->NHits = hits->size();

  int minLayer(34);
  int maxLayer(-1);

  // Loop through the hits.
  for(COMET::IHitSelection::iterator hit = hits->begin(); hit != hits->end(); ++hit){

    int layer = COMET::IGeomInfo::Get().ECAL().GetLayerNumber((*hit)->GetGeomId());

    if(layer < minLayer) minLayer = layer;
    if(layer > maxLayer) maxLayer = layer;

    this->TotalHitCharge += (*hit)->GetCharge();
    this->AverageHitTime += (*hit)->GetTime() * (*hit)->GetCharge();
  }

  // Finish charge weighting the time.
  this->AverageHitTime /= this->TotalHitCharge;

  // Layer Info
  this->mostUpStreamLayerHit = minLayer;
  this->mostDownStreamLayerHit = maxLayer;

  double minLayerCharge(0.0);
  double maxLayerCharge(0.0);

  // Find Charge Weighted Front and Back Positions
  for(COMET::IHitSelection::iterator hit = hits->begin(); hit != hits->end(); ++hit){
    TVector3 hitPos = (*hit)->GetPosition();
    int layer = COMET::IGeomInfo::Get().ECAL().GetLayerNumber((*hit)->GetGeomId());

    if(layer == minLayer){
      this->FrontPosition += hitPos * (*hit)->GetCharge();
      minLayerCharge += (*hit)->GetCharge();
    }
    if(layer == maxLayer){
      this->BackPosition += hitPos * (*hit)->GetCharge();
      maxLayerCharge += (*hit)->GetCharge();
    }
  }
  this->FrontPosition *= (1.0 / minLayerCharge);
  this->BackPosition *= (1.0 / maxLayerCharge);

  // Find the view from one of the layers used.
  if(minLayer%2 == 0) this->View = 1;
  else this->View = 0;

    this->G4ID = TrackTruthInfo::GetG4TrajIDHits(*hits);
}

void COMET::ITrackerReconECalModule::IReconECalObject::CalculateFrontBackRatio(COMET::IHandle<COMET::IHitSelection> hits, double& FrontBackRatio, double& FrontBackRatioNumerator, double& FrontBackRatioDenominator)
{
  //create a vector to store the hits in
  std::vector< std::pair<TVector3,double> > positionAndCharge;

  if (hits && hits->size())
  {
    //determine principal PCA axis
    TPrincipal pca(3,"");

    for(COMET::IHitSelection::iterator it = hits->begin(); it != hits->end(); ++it)
    {
      const TVector3& pos = (*it)->GetPosition();
      double charge = (*it)->GetCharge();
      int repeatEntries = static_cast<int>(charge / 0.1);
      for(int i = 0; i < repeatEntries; ++i)
      {
        double row[3] = { pos.X(), pos.Y(), pos.Z() };
        pca.AddRow(row);
      }
    }
    pca.MakePrincipals();
    
    //Rotate hits into PCA space, and fill vector with the new positions
    for(COMET::IHitSelection::iterator it = hits->begin(); it != hits->end(); ++it)
    {
      TVector3 position = (*it)->GetPosition();
      COMET::IHandle<COMET::IHit> h = *it;
        double charge = h->GetCharge();
      
      double x[3] = { position.X(),
              position.Y(),
              position.Z()};
      double p[3];
      pca.X2P(x,p);
      position.SetZ(p[0]); //define Z-axis of new co-ordinate system to be the direction along the principal PCA axis
      position.SetY(p[1]);
      position.SetX(p[2]);

      std::pair<TVector3, double> p_c(position, charge);
        positionAndCharge.push_back( p_c );
    }

    //define four blocks, first to fourth
    double first(0);
    double second(0);
    double third(0);
    double fourth(0);
    
    
    int size = positionAndCharge.size();
    //set initial values for the ends
    double zmin = positionAndCharge.at(0).first.Z();
    double zmax = positionAndCharge.at(0).first.Z();
    
    //find the actual maximum and minimum "z" values (i.e. the ends of the track)
    for (int i(0); i < size ; i++)
    {
      if (positionAndCharge.at(i).first.Z() > zmax)
      {
        zmax = positionAndCharge.at(i).first.Z();
      }
    
      if (positionAndCharge.at(i).first.Z() < zmin)
      {
        zmin = positionAndCharge.at(i).first.Z();
      }
    }

    //calculate the lengh of the track along the PCA axis
    double length = zmax - zmin;

    //loop through the hits, calculating which block it falls in, and adding its charge to that block
    for(int i = 0; i < size; i++)
    {
        if ( (positionAndCharge.at(i).first.Z() - zmin) > (3 * length/4) )
      {
            fourth += positionAndCharge.at(i).second;
      }
      else if ( (positionAndCharge.at(i).first.z() - zmin) > (length/2))
      {
        third += positionAndCharge.at(i).second;
      }
      else if ( (positionAndCharge.at(i).first.z() - zmin) > (length/4) )
      {
        second += positionAndCharge.at(i).second;
      }
      else
      {
        first += positionAndCharge.at(i).second;
      }
    }
    
    //if the filling has succeeded, calculate the ratio
    if (first != 0 && fourth != 0)
    {
      FrontBackRatioNumerator = fourth;
      FrontBackRatioDenominator = first;
      FrontBackRatio = fourth / first;
    }
  }
}

COMET::ITrackerReconECalModule::IReconECalTrack::~IReconECalTrack(){}
COMET::ITrackerReconECalModule::IReconECalShower::~IReconECalShower(){}
COMET::ITrackerReconECalModule::IReconECalCluster::~IReconECalCluster(){}
COMET::ITrackerReconECalModule::IReconECalObject::~IReconECalObject(){}
COMET::ITrackerReconECalModule::IReconECalUnmatchedObject::~IReconECalUnmatchedObject(){}