MakeNoiseSpectrumFile_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   TODO
/// \author  Christopher Backhouse
////////////////////////////////////////////////////////////////////////

#include "ReadoutSim/FPGAAlgorithms.h"
#include "ReadoutSim/NoiseMaker.h"

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"

#include "fhiclcpp/ParameterSet.h"

#include "CLHEP/Random/RandGaussQ.h"

// NOvA includes
#include "Simulation/Simulation.h"

#include "TH3.h"
#include "TMatrixD.h"
#include "TVector.h"

namespace rsim
{
  const int kCovMatSize = 6;

  class MakeNoiseSpectrumFile: public art::EDAnalyzer, protected CommonParameters
  {
  public:
    explicit MakeNoiseSpectrumFile(fhicl::ParameterSet const &pset);
    virtual ~MakeNoiseSpectrumFile();

    void analyze(const art::Event& evt);
    void beginJob();
    void endJob();

  protected:
    int fN; ///< How many spills have we simulated?

    // Information for calculating the covariance matrix of the noise.
    // Need to use TMatrixD because default TMatrix overflows.
    TMatrixD fCovMat;
    TMatrixD fCovMatStats;
    TVectorD fMeanVec;
    TVectorD fMeanVecStats;

    rsim::IFPGAAlgorithm* fFPGAAlgo;
    rsim::NoiseMaker    fNoiseMaker;

    TH3* fNoiseMap;
  };

  //----------------------------------------------------------------------
  MakeNoiseSpectrumFile::MakeNoiseSpectrumFile(fhicl::ParameterSet const& pset) :
    EDAnalyzer(pset),
    CommonParameters(pset),
    fN(0),
    fCovMat(kCovMatSize, kCovMatSize), fCovMatStats(kCovMatSize, kCovMatSize),
    fMeanVec(kCovMatSize), fMeanVecStats(kCovMatSize),
    fFPGAAlgo(CreateFPGAAlgorithm(pset)),
    fNoiseMaker(pset)
  {
    // get the random number seed, use a random default if not specified
    // in the configuration file.
    unsigned int seed = pset.get<unsigned int>("Seed",
                                               sim::GetRandomNumberSeed());

    createEngine(seed);
  }

  //----------------------------------------------------------------------
  MakeNoiseSpectrumFile::~MakeNoiseSpectrumFile()
  {
    delete fFPGAAlgo;
  }

  //----------------------------------------------------------------------
  void MakeNoiseSpectrumFile::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;

    fNoiseMap = tfs->make<TH3D>("h3_adc0adc1adc2", ";ADC0;ADC1;ADC2",
                                80, 0, 80,
                                80, 0, 80,
                                80, 0, 80);
  }

  //----------------------------------------------------------------------
  void MakeNoiseSpectrumFile::analyze(const art::Event& evt)
  {
    ++fN;

    art::ServiceHandle<art::RandomNumberGenerator> rng;
    std::cout.flush();
    CLHEP::RandGaussQ gauss(rng->getEngine());

    const int maxADC = (1<<12)-1;
    double ASIC_Output[fNumClockticksInSpill];

    for(int n = 0; n < fNumClockticksInSpill;
        n += fNumClockticksPerDigitization){
      ASIC_Output[n] = int(fASICBaselineInADCCounts +
                           fNoiseMaker.GetNoise(&gauss)/fADCMaxPE*maxADC);

      if(ASIC_Output[n] < 0) ASIC_Output[n] = 0;
      if(ASIC_Output[n] > maxADC) ASIC_Output[n] = maxADC;
    }


    // Fill covariance matrix information
    for(int n = 0; n < fNumClockticksInSpill; n += fNumClockticksPerDigitization){
      for(int i = 0; i < kCovMatSize; ++i){
        const int iidx = n+fNumClockticksPerDigitization*i;
        if(iidx >= fNumClockticksInSpill) continue;

        fMeanVec(i) += ASIC_Output[iidx];
        ++fMeanVecStats(i);

        for(int j = 0; j < kCovMatSize; ++j){
          const int jidx = n+fNumClockticksPerDigitization*j;
          if(jidx >= fNumClockticksInSpill) continue;

          fCovMat(i, j) += ASIC_Output[iidx]*ASIC_Output[jidx];
          ++fCovMatStats(i, j);
        } // end for j
      } // end for i
    } // end for n


    const std::vector<rawdata::RawDigit> digits =
      fFPGAAlgo->ASICToDigits(ASIC_Output, 0);

    for(unsigned int i = 0; i < digits.size(); ++i){
      const rawdata::RawDigit& digit = digits[i];

      fNoiseMap->Fill(digit.ADC(0), digit.ADC(1), digit.ADC(2));
    }
  }

  //----------------------------------------------------------------------
  void MakeNoiseSpectrumFile::endJob()
  {
    // Normalization of noise map is probability /cell /ns
    // Warning, this can overflow if performed with ints.
    const double ns = long(fN)*long(fNumClockticksInSpill)*fClocktick;
    fNoiseMap->Scale(1/ns);

    // Divide by number of samples to get means
    for(int i = 0; i < kCovMatSize; ++i){
      fMeanVec(i) /= fMeanVecStats(i);
      for(int j = 0; j < kCovMatSize; ++j){
        fCovMat(i, j) /= fCovMatStats(i, j);
      }
    }

    // Definition of covariance: <ij>-<i><j>
    for(int i = 0; i < kCovMatSize; ++i){
      for(int j = 0; j < kCovMatSize; ++j){
        fCovMat(i, j) -= fMeanVec(i)*fMeanVec(j);
      }
    }

    std::cout << "Covariance matrix: " << std::endl;
    fCovMat.Print();
    std::cout << "And its inverse: " << std::endl;
    fCovMat.Invert();
    fCovMat.Print();

    /*
    // Output format useful for copy-pasting into code
    for(int i = 0; i < kCovMatSize; ++i){
      for(int j = 0; j < kCovMatSize; ++j){
        std::cout << "mat(" << i << ", " << j << ") = "
                  << fCovMat(i, j) << "; ";
      }
      std::cout << std::endl;
    }
    */
  }

  DEFINE_ART_MODULE(MakeNoiseSpectrumFile)

} // namespace