////////////////////////////////////////////////////////////////////
/// \class RAT::DS::MC
///
/// \brief  This class represents a single simulated Monte Carlo event.
///
/// \author Stan Seibert <volsung@physics.utexas.edu>\n
///         Phil G Jones <p.g.jones@qmul.ac.uk>
///
/// REVISION HISTORY:\n
///   8 July 2009 : W. Heintzelman - Expand track pruning capability so
///                 tracks can be "condensed"\n
///  22 Feb  2010 : W. Heintzelman - Expand track pruning capability so tracks
///                             can be pruned based on the type of
///                             physical process that produced them.\n
///  02 Apr 2010 : Gabriel Orebi Gann - Add real time to MC (i.e. date etc)\n
///  28 Sep 2012 : A. Mastbaum - Remove Braidwood muon veto DS\n
///  2013-1-21 : P. Jones - Refactor as part of ds review.
///  2016-10-05  : W Heintzelman - Add  GetMCTrackFromIndex and  GetMCTrackID
///                                functions; correct comments
///                                Add numDirPE variable and set/get functions
///  2016-10-19 : M. Stringer-MC Tracks now indexed by their trackID. (PR #1508)
///
/// \details This structure contains the "Monte Carlo Truth" for the event.
/// Intermediate simulated quantities, like electronics response, can
/// also go here.  Only once the data finally looks like "real life"
/// should it go into the EV class.
///
////////////////////////////////////////////////////////////////////
#ifndef __RAT_DS_MC__
#define __RAT_DS_MC__

#include <TObject.h>

#include <RAT/ListHelp.hh>
#include <RAT/DS/UniversalTime.hh>
#include <RAT/DS/MCParticle.hh>
#include <RAT/DS/MCTrack.hh>
#include <RAT/DS/MCPMT.hh>
#include <RAT/DS/MCHit.hh>
#include <RAT/DS/PMTSet.hh>

#include <vector>
#include <set>
#include <algorithm>

namespace RAT
{
namespace DS
{

class MC : public TObject
{
public:
  MC() : TObject(), mcTime(0), mcid(0), mcPECount(0), nCherPhotons(0), nScintPhotons(0),
         nReemittedPhotons(0), intialScintTime(0), scintEnergyDeposit(0), scintQuenchedEnergyDeposit(0) { };

  /// Set the event number/ID
  ///
  /// @param[in] eventID for the current event
  void SetMCID( const Int_t eventID ) { mcid = eventID; }

  /// Get the event number/ID
  ///
  /// @return the event ID
  Int_t GetMCID() const { return mcid; }

  /// Set the event Monte Carlo, MC, time
  ///
  /// This is the time of the event in ns since the beginning of the simulation
  ///
  /// @param[in] time of the event relative the MC start
  void SetMCTime( const Double_t time ) { mcTime = time; }

  /// Return the event Monte Carlo, MC, time
  ///
  /// This is the time of the event in ns since the beginning of the simulation
  ///
  /// @return time of the event relative the MC start
  Double_t GetMCTime() const { return mcTime; }

  /// Set the event Universal Time
  ///
  /// This is the time of the event relative to the SNO+ time origin
  ///
  /// @param[in] time of the event in the universal time system
  void SetUniversalTime( const UniversalTime& time ) { universalTime = time; }

  /// Get the event Universal Time
  ///
  /// This is the time of the event relative to the SNO+ time origin
  ///
  /// @return time of the event in the universal time system
  UniversalTime GetUniversalTime() const { return universalTime; }

  /// Add a new parent particle
  ///
  /// @param[in] particle the parent particle information
  void AddMCParent( const MCParticle& particle ) { parents.push_back( particle ); }

  /// Get the total number of parent particles
  ///
  /// @return the number of parent particles
  size_t GetMCParentCount() const { return parents.size(); }

  /// Get a MCParent using the index in the parents vector
  ///
  /// @param[in] index into the parents array
  /// @return reference to the MCParent
  /// @throws out_of_range if the index is out of range
  MCParticle& GetMCParent( const size_t index ) { return parents.at( index ); }
  /// @copydoc GetMCParent(size_t)
  const MCParticle& GetMCParent( const size_t index ) const { return parents.at( index ); }

  /// Prune/delete the MCParent collection
  void PruneMCParents() { clear_vector<MCParticle>( parents ); }

  /// Add a new particle
  ///
  /// @param[in] particle the particle information
  void AddMCParticle( const MCParticle& particle ) { particles.push_back( particle ); }

  /// Get the total number of particles
  ///
  /// @return the number of particles
  size_t GetMCParticleCount() const { return particles.size(); }

  /// Get a MCParticle using the index in the particles vector
  ///
  /// @param[in] index into the particles array
  /// @return reference to the MCParticle
  /// @throws out_of_range if the index is out of range
  MCParticle& GetMCParticle( const size_t index ) { return particles.at( index ); }
  /// @copydoc GetMCParticle(size_t)
  const MCParticle& GetMCParticle( const size_t index ) const { return particles.at( index ); }

  /// Prune/delete the MCParticle collection
  void PruneMCParticles() { clear_vector<MCParticle>( particles ); }

  /// Add a new track
  ///
  /// @param[in] track the track information
  void AddMCTrack( const MCTrack& track ) {
      tracks.push_back( track );
      trackIDs.push_back(track.GetTrackID());
  }

  /// Get the total number of tracks
  ///
  /// Tracks only exist if
  /// @verbatim
  /// /tracking/storeTrajectory 1
  /// @endverbatim
  /// is present in the macro
  ///
  /// @return the number of tracks
  size_t GetMCTrackCount() const { return tracks.size(); }

  /// Get a vector of track IDs
  ///
  /// Tracks only exist if
  /// @verbatim
  /// /tracking/storeTrajectory 1
  /// @endverbatim
  /// is present in the macro
  ///
  /// @return a vector containing all the trackIDs
  std::vector<unsigned int> GetMCTrackIDs() const { return trackIDs; }

  /// Get a MCTrack using the track's ID
  ///
  /// Tracks only exist if
  /// @verbatim
  /// /tracking/storeTrajectory 1
  /// @endverbatim
  /// is present in the macro
  ///
  /// @param[in] ID of the track
  /// @return reference to the MCTrack
  /// @throws out_of_range if the index is out of range
  MCTrack& GetMCTrack( const size_t trackID ){
      for(unsigned i=0; i<tracks.size(); i++){
          if(trackIDs.at(i) == trackID){
              return tracks.at(i);
          }
      }
      //If we can't find the track try and return somthing out of bounds to throw out_of_range exception
      return tracks.at(tracks.size()+1);
  }

  /// @copydoc GetMCTrack(size_t)
  const MCTrack& GetMCTrack( const size_t trackID ) const {
      for(unsigned i=0; i<tracks.size(); i++){
          if(trackIDs.at(i) == trackID){
              return tracks.at(i);
          }
      }
      //If we can't find the track try and return somthing out of bounds to throw out_of_range exception
      return tracks.at(tracks.size()+1);
  }

  /// Get a MCTrack using the index in the tracks vector
  MCTrack& GetMCTrackFromIndex(const size_t index ){ return tracks.at( index );}

  /// Get the MC Track ID given the index in the tracks vector
  size_t GetMCTrackID( const size_t index ){
      return tracks.at(index).GetTrackID();
  }

  /// Prune/delete the MCTrack collection
  void PruneMCTracks() { clear_vector<MCTrack>( tracks ); }

  /// Prune specific tracks by index in the MCTrack collection
  ///
  /// @param[in] indices of the tracks to prune
  inline void PruneMCTracksByIndex( const std::set<size_t>& indices );

  /// Add a new pmt
  ///
  /// @param[in] pmt the pmt information
  void AddMCPMT( const MCPMT& pmt ) { pmts.push_back( pmt ); }

  /// Get the total number of pmts
  ///
  /// @return the number of pmts
  size_t GetMCPMTCount() const { return pmts.size(); }

  /// Get a MCPMT using the index in the pmts vector
  ///
  /// @param[in] index into the pmts array
  /// @return reference to the MCPMT
  /// @throws out_of_range if the index is out of range
  MCPMT& GetMCPMT( const size_t index ) { return pmts.at( index ); }
  /// @copydoc GetMCPMT(size_t)
  const MCPMT& GetMCPMT( const size_t index ) const { return pmts.at( index ); }

  /// Prune/delete the MCPMT collection
  void PruneMCPMTs() { clear_vector<MCPMT>( pmts ); }

  /// Get the Unbuilt MCHit collection
  ///
  /// These are Hits that are above the channel threshold but are not associated
  /// with or built into a detector event (MCEV)
  ///
  /// @return reference to the MCHits
  MCHits& GetUnbuiltMCHits() { return unbuiltHits; }
  /// @copydoc GetUnbuiltMCHits()
  const MCHits& GetUnbuiltMCHits() const { return unbuiltHits; }

  /// Set the number of photoelectrons in this event
  ///
  /// @param[in] photoelectrons total number of
  void SetMCPECount( const UInt_t photoelectrons ) { mcPECount = photoelectrons; }

  /// Get the number of photoelectrons in this event
  ///
  /// @return number of photoelectrons
  UInt_t GetMCPECount() const { return mcPECount; }

  /// Total number of photoelectrons generated by photons whose origin
  /// is within 1.1*(concentrator radius) of the PMT position in which
  /// the photoelectron was generated  (calculated in EventBuilderProc.cc)
  virtual void SetNumDirPE(UInt_t _numDirPE) { numDirPE = _numDirPE; };
  virtual UInt_t GetNumDirPE() const { return numDirPE; };

  /// Set the number of Cerenkov Photons in this event
  ///
  /// @param[in] cerenkovPhotons total number of
  void SetNCherPhotons( const UInt_t cerenkovPhotons ) { nCherPhotons = cerenkovPhotons; }

  /// Get the number of Cerenkov Photons in this event
  ///
  /// @return number of Cerenkov Photons
  UInt_t GetNCherPhotons() const { return nCherPhotons; }

  /// Set the number of Scintillated Photons in this event
  ///
  /// @param[in] scintillatedPhotons total number of scintillated
  void SetNScintPhotons( const UInt_t scintillatedPhotons ) { nScintPhotons = scintillatedPhotons ; }

  /// Get the number of Scintillated Photons in this event
  ///
  /// @return number of Scintillated Photons
  UInt_t GetNScintPhotons() const { return nScintPhotons; }

  /// Set the number of reemitted Photons in this event
  ///
  /// @param[in] reemittedPhotons total number of reemitted photons
  void SetNReemittedPhotons( const UInt_t reemittedPhotons ) { nReemittedPhotons = reemittedPhotons; }

  /// Get the number of reemitted Photons in this event
  ///
  /// @return number of reemitted Photons
  UInt_t GetNReemittedPhotons() const { return nReemittedPhotons; }

  /// Set the initial scintillation time
  ///
  /// This is the initial time of scintillation
  ///
  /// @param[in] initialTime of scintillation
  void SetInitialScintTime( const Double_t initialTime ) { intialScintTime = initialTime; }

  /// Get the initial scintillation time
  ///
  /// This is the initial time of scintillation
  ///
  /// @return initial time of scintillation
  Double_t GetInitialScintTime() const { return intialScintTime; }

  /// Set the total energy deposited in the scintillator
  ///
  /// This is energy deposited in any active scintillator
  ///
  /// @param[in] energy deposited in the scintillator
  void SetScintEnergyDeposit( const Double_t energy ) { scintEnergyDeposit = energy; }

  /// Get the total energy deposited in the scintillator
  ///
  /// This is energy deposited in any active scintillator
  ///
  /// @return energy deposited in the scintillator
  Double_t GetScintEnergyDeposit() const { return scintEnergyDeposit; }

  /// Set the total quenched energy deposited in the scintillator
  ///
  /// @param[in] energy quenched deposited in the scintillator
  void SetScintQuenchedEnergyDeposit( const Double_t energy ) { scintQuenchedEnergyDeposit = energy; }

  /// Get the total quenched energy deposited in the scintillator
  ///
  /// @return quenched energy deposited in the scintillator
  Double_t GetScintQuenchedEnergyDeposit() const { return scintQuenchedEnergyDeposit; }

  // This ROOT macro adds dictionary methods to this class.
  // The number should be incremented whenever this class's members are changed.
  // It assumes this class has no virtual methods, use ClassDef if change this.
  ClassDefNV( MC, 11 );
protected:
  std::vector<MCParticle> particles; ///< The primary particles for this event
  std::vector<MCParticle> parents; ///< Primary particle parent particles e.g. neutrino or decay parent
  std::vector<MCTrack> tracks; ///< Particle tracks as simulated (empty by default)
  std::vector<unsigned int> trackIDs; ///< Vector of Track IDs(empty by default)
  std::vector<MCPMT> pmts; ///< PMTs that record at least one photoelectron in the simulation before the daq, or a photon in the concentrator bucket
  MCHits unbuiltHits; ///< Unbuilt MCHits, hits not associated with a detector event

  UniversalTime universalTime; ///< Universal time of the event
  Double_t mcTime; ///< MC time i.e. relative to beginning of run in ns

  Int_t mcid; ///< The ID of this event

  UInt_t mcPECount; ///< Total number of photoelectrons (before the daq)
  UInt_t nCherPhotons; ///< Number of cerenkov photons created in the event
  UInt_t nScintPhotons; ///< Number of scintillation photons created in the event
  UInt_t nReemittedPhotons; ///< Number of reemitted photons in the event
  UInt_t numDirPE;          ///< Used with GenPMTEff event generator

  Double_t intialScintTime; ///< Initial time of scintillation (t_zero)
  Double32_t scintEnergyDeposit; ///< Total energy deposited in the scintillator
  Double32_t scintQuenchedEnergyDeposit; ///< Total visible/quenched energy deposited in the scintillator
};

inline void
MC::PruneMCTracksByIndex( const std::set<size_t>& indices )
{
    for( std::set<size_t>::reverse_iterator iTer = indices.rbegin(); iTer != indices.rend(); iTer++ ){
        #ifdef RATVERBOSE
          detail << "[MC]::PruneMCTracksByIndex : Cutting track: "<<*iTer<<" which is an: "<<tracks.at(*iTer).GetParticleName()<<newline;
        #endif
        tracks.erase(tracks.begin()+*iTer);
        trackIDs.erase(trackIDs.begin()+*iTer);
      }
}

} // namespace DS

} // namespace RAT

#endif