#include <globals.hh>
#include <Randomize.hh>

#include <RAT/Noise.hh>
#include <RAT/Log.hh>
#include <RAT/DB.hh>
#include <RAT/ChannelEfficiency.hh>
#include <RAT/HitPMTCollection.hh>
#include <RAT/DS/MCPMT.hh>
#include <RAT/DS/MCPE.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/ChanHWStatus.hh>
#include <RAT/DU/PMTInfo.hh>
using namespace RAT;

#include <CLHEP/Units/SystemOfUnits.h>
#include <vector>
using namespace CLHEP;
using namespace std;

Noise::Noise()
{

}

Noise::~Noise()
{

}

void
Noise::BeginOfRun()
{
  DBLinkPtr daqDB = DB::Get()->GetLink( "DAQ" );
  fTriggerWindow = daqDB->GetD( "triggergate" );

  DBLinkPtr noiseDB = DB::Get()->GetLink( "NOISE_RUN_LEVEL" );
  DBLinkPtr integrated_noiseDB = DB::Get()->GetLink( "NOISE_RUN_INTEGRATED" );
  DBLinkPtr noiseMC = DB::Get()->GetLink( "NOISE_MC");

  fNoiseFlag = static_cast<ENoiseModel>( noiseMC->GetI( "noise_flag" ) );
  fTableFlag = static_cast<ENoiseTable>( noiseMC->GetI( "integ_noise_flag" ) );

  if( fTableFlag == eRunLevelTable )
    fNoiseRate = noiseDB->GetD( "average_noiserate" ) / s;     // Rate per ns
  else if( fTableFlag == eIntegratedTable )
    fNoiseRate = integrated_noiseDB->GetD( "average_noiserate" ) / s;
  else
    Log::Die( "Noise::BeginOfRun Table flag is not valid" );

  fNoiseGroupFlag = noiseMC->GetI( "noise_group_flag" ) == 1;
  fNoiseTEarly = noiseMC->GetD( "noise_group_early" );
  fNoiseTLate = noiseMC->GetD( "noise_group_late" );

  if( fNoiseFlag == eNoNoise )
    info << "Noise::BeginOfRun PMT noise: off" << newline;
  else if( fNoiseFlag == eAverageNoise )
    info << "Noise::BeginOfRun PMT noise: average rate" << newline;
  else if( fNoiseFlag == ePerChannelNoise )
    info << "Noise::BeginOfRun PMT noise: per-PMT rates from db" << newline;
  else
    Log::Die( "Noise::BeginOfRun Noise flag is not valid" );
  if( fNoiseGroupFlag )
    info << "Noise::BeginOfRun grouping enabled." << newline;

  // Get tube status and noise rate
  fNumChannels = DB::Get()->GetLink( "PMTINFO" )->GetIArray( "pmt_type" ).size();
  if( fTableFlag == eRunLevelTable )
    fChannelRate = noiseDB->GetDArray( "perpmt_noiserate" );
  else if( fTableFlag == eIntegratedTable )
    fChannelRate = integrated_noiseDB->GetDArray( "perpmt_noiserate" );
  else
    Log::Die( "Noise::BeginOfRun Table flag is not valid" );
  for( vector<double>::iterator iChannelRate = fChannelRate.begin(); iChannelRate != fChannelRate.end(); ++iChannelRate )
    *iChannelRate /= s; // Rate in seconds

  fDefaultWindowStart = -200.0 * ns;
  fDefaultWindowEnd = 600.0 * ns;

  fOWLnoise = noiseMC->GetI("OWLnoise");
  fAveChanEff = 0.;
  if(fOWLnoise==0) fAveChanEff = ChannelEfficiency::Get()->GetAverageEfficiency();
  if(fOWLnoise==1) fAveChanEff = ChannelEfficiency::Get()->GetAverageEffPlusOWLs();


}

void
Noise::GenerateNoise()
{
  if( fNoiseFlag == eNoNoise )
    // Drop out, no noise to generate
    return;

  // Get the collection of Monte Carlo hits
  const std::vector<DS::MCPMT*>& pmtHits = HitPMTCollection::Get()->GetSortedPMTCollection();
  // Firstly generate noise in the absence of hits over the full default window
  if( pmtHits.size() == 0 )
    {
      GenerateNoiseInWindow( fDefaultWindowStart, fDefaultWindowEnd );
      // Drop out, noise generated
      return;
    }
  // Add noise alongside existing hits, first need to know an ordered list of hit times.
  // This is required to know the first and last hit times and allow for grouping of
  // hits.
  vector<double> hitTimes;
  for( std::vector<DS::MCPMT*>::const_iterator iTer = pmtHits.begin(); iTer != pmtHits.end(); ++iTer )
    {
      for( size_t iPhotoelectron = 0; iPhotoelectron < (*iTer)->GetMCPECount(); iPhotoelectron++ )
        hitTimes.push_back( (*iTer)->GetMCPE( iPhotoelectron ).GetCreationTime() );
    }
  // Sort the hit times
  std::sort( hitTimes.begin(), hitTimes.end() );
  // Can now add noise, firstly see if we need to group the hits, if not just generate
  const double windowStart = hitTimes[0];
  const double windowEnd = hitTimes[hitTimes.size() - 1];
  if( fNoiseGroupFlag == false )
    {
      GenerateNoiseInWindow( windowStart, windowEnd );
      // Drop out, noise generated
      return;
    }
  // Must group the hits
  else
    {
      double groupStart = windowStart;
      double groupEnd = windowEnd;
      for( size_t iHit = 1; iHit < hitTimes.size(); iHit++ )
        {
          // If spacing to previous hit is larger than the grouping window, then group and generate
          if( hitTimes[iHit] > hitTimes[iHit - 1] + fTriggerWindow + fNoiseTEarly + fNoiseTLate )
            {
              groupEnd = hitTimes[iHit - 1];
              GenerateNoiseInWindow( groupStart, groupEnd );
              groupStart = hitTimes[iHit];
            }
        }
      // Generate noise in the final group
      GenerateNoiseInWindow( groupStart, windowEnd );
    }
}

void
Noise::GenerateNoiseInWindow( double windowStart, ///< Start MC time in ns
                              double windowEnd ) ///< End MC time in ns
{
  windowStart -= fNoiseTEarly; // Early buffer
  windowEnd += fNoiseTLate + fTriggerWindow; // Late buffer (includes trigger window length)
  switch( fNoiseFlag )
    {
    case eAverageNoise:
      GenerateAverageNoise( windowStart, windowEnd );
      break;
    case ePerChannelNoise:
      GeneratePerChannelNoise( windowStart, windowEnd );
      break;
    default:
      break;
    }
}

void
Noise::GenerateAverageNoise( const double windowStart,
                             const double windowEnd )
{
  // Calculate the average noise rate per detector corrected by the AVERAGE detection efficiency
  // First calculate the expected number of noise hits, then generate these hits on random channels
  // NOTE: the expected number of hits, and the channels chosen, are at this point across the whole detector
  // (this is self consistent)
  // i.e. not worrying about online channels or normal channels
  // This will be handled in the GenerateNoiseHit funtion, which will only generate the hit for online, normal/HQE/OWL channels

  const double perChannelNoise = fNoiseRate * ( windowEnd - windowStart );
  const double perDetectorNoise = perChannelNoise * fNumChannels / fAveChanEff;
  const int numNoiseHits = static_cast<int>( floor( CLHEP::RandPoisson::shoot( perDetectorNoise ) ) );
  for( int iNoise = 0; iNoise < numNoiseHits; iNoise++ )
    {
      int lcn = static_cast<int>( G4UniformRand() * fNumChannels );
      GenerateNoiseHit( lcn, windowStart, windowEnd );
    }
}

void
Noise::GeneratePerChannelNoise( const double windowStart,
                                const double windowEnd )
{
  // Calculate the average noise rate for this channel corrected by the CHANNEL detection efficiency
  // First calculate the expected number of noise hits for each channel individually, then generate
  // noise hits.
  for( int lcn = 0; lcn < fNumChannels; lcn++ )
    {
      double channelEfficiency = ChannelEfficiency::Get()->GetChannelEfficiency( lcn );
      // If no channel efficiency is known, set it to the average
      if( channelEfficiency == 0.0 )
        channelEfficiency = fAveChanEff;
      const double noiseRate = fChannelRate[lcn] * ( windowEnd - windowStart ) / channelEfficiency;
      const int numNoiseHits = static_cast<int>( floor( CLHEP::RandPoisson::shoot( noiseRate ) ) );
      for( int iNoise = 0; iNoise < numNoiseHits; iNoise++ )
        GenerateNoiseHit( lcn, windowStart, windowEnd );
    }
}

void
Noise::GenerateNoiseHit( const int lcn,
                         const double windowStart,
                         const double windowEnd )
{
  const DU::ChanHWStatus& channelHardwareStatus = DU::Utility::Get()->GetChanHWStatus();
  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();

  // Only generate noise on online channels (if we're checking DQXX)
  if( channelHardwareStatus.IsEnabled() && !channelHardwareStatus.IsTubeOnline( lcn ) ) return;

  // Only generate noise on normal and HQE tubes, and OWLs if selected
  if( !pmtInfo.GetType(lcn)==DU::PMTInfo::NORMAL && !pmtInfo.GetType(lcn)==DU::PMTInfo::HQE && (fOWLnoise==0 || !pmtInfo.GetType(lcn)==DU::PMTInfo::OWL )) return;

  // Choose a random time in this window
  const double time = windowStart + G4UniformRand() * ( windowEnd - windowStart );
  // create new DS::MCPE, the way of recording photo hits on PMTs
  DS::MCPE photoelectron;
  photoelectron.SetCreationTime( time );
  photoelectron.SetCharge( 1 );
  photoelectron.SetNoise( true );
  HitPMTCollection::Get()->DetectPhotoelectron( lcn, photoelectron );
}