IDCFourierTransform.h

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#ifndef ALICEO2_IDCFOURIERTRANSFORM_H_
#define ALICEO2_IDCFOURIERTRANSFORM_H_
 
#include <vector>
#include "Rtypes.h"
#include "DataFormatsTPC/Defs.h"
#include "TPCCalibration/IDCContainer.h"
#include "TPCCalibration/IDCFourierTransformBase.h"
#include "CommonConstants/LHCConstants.h"
 
using fftwf_plan = struct fftwf_plan_s*;
using fftwf_complex = float[2];
 
namespace o2::tpc
{
 
 
template <class Type> // do not use enum class as type to avoid problems with ROOT dictionary generation!
class IDCFourierTransform : public IDCFourierTransformBase<Type>
{
 public:
  template <bool IsEnabled = true, typename std::enable_if<(IsEnabled && (std::is_same<Type, IDCFourierTransformBaseAggregator>::value)), int>::type = 0>
  IDCFourierTransform(const unsigned int rangeIDC = 200, const unsigned int nFourierCoefficientsStore = 200 + 2) : IDCFourierTransformAggregator(rangeIDC), mFourierCoefficients{1, nFourierCoefficientsStore}, mVal1DIDCs(sNThreads), mCoefficients(sNThreads)
  {
  }
 
  template <bool IsEnabled = true, typename std::enable_if<(IsEnabled && (std::is_same<Type, IDCFourierTransformBaseEPN>::value)), int>::type = 0>
  IDCFourierTransform(const unsigned int rangeIDC = 200, const unsigned int nFourierCoefficientsStore = 200 + 2) : IDCFourierTransformEPN(rangeIDC), mFourierCoefficients{1, nFourierCoefficientsStore}, mVal1DIDCs(sNThreads), mCoefficients(sNThreads)
  {
  }
 
  // Destructor
  ~IDCFourierTransform();
 
  static void setFFT(const bool fft) { sFftw = fft; }
 
  template <bool IsEnabled = true, typename std::enable_if<(IsEnabled && (std::is_same<Type, IDCFourierTransformBaseAggregator>::value)), int>::type = 0>
  static void setNThreads(const int nThreads)
  {
    sNThreads = nThreads;
  }
 
  void initFFTW3Members();
 
  template <bool IsEnabled = true, typename std::enable_if<(IsEnabled && (std::is_same<Type, IDCFourierTransformBaseAggregator>::value)), int>::type = 0>
  void calcFourierCoefficients(const unsigned int timeFrames = 2000)
  {
    mFourierCoefficients.resize(timeFrames);
    sFftw ? calcFourierCoefficientsFFTW3() : calcFourierCoefficientsNaive();
  }
 
  template <bool IsEnabled = true, typename std::enable_if<(IsEnabled && (std::is_same<Type, IDCFourierTransformBaseEPN>::value)), int>::type = 0>
  void calcFourierCoefficients()
  {
    sFftw ? calcFourierCoefficientsFFTW3() : calcFourierCoefficientsNaive();
  }
 
  std::vector<std::vector<float>> inverseFourierTransform() const { return sFftw ? inverseFourierTransformFFTW3() : inverseFourierTransformNaive(); }
 
  unsigned int getrangeIDC() const { return this->mRangeIDC; }
 
  const auto& getFourierCoefficients() const { return mFourierCoefficients; }
 
  static bool getFFT() { return sFftw; }
 
  static int getNThreads() { return sNThreads; }
 
  void dumpToFile(const char* outFileName = "Fourier.root", const char* outName = "FourierCoefficients") const;
 
  void dumpToTree(const char* outFileName = "FourierTree.root") const;
 
  void printFFTWPlan() const;
 
  std::vector<std::pair<float, float>> getFrequencies(const float samplingFrequency = getSamplingFrequencyIDCHz()) const { return getFrequencies(mFourierCoefficients, samplingFrequency); }
 
  static std::vector<std::pair<float, float>> getFrequencies(const FourierCoeff& coeff, const float samplingFrequency = getSamplingFrequencyIDCHz());
 
  static float getSamplingFrequencyIDCHz() { return 1e6 / (12 * o2::constants::lhc::LHCOrbitMUS); }
 
 private:
  FourierCoeff mFourierCoefficients;         
  inline static int sFftw{1};                
  inline static int sNThreads{1};            
  fftwf_plan mFFTWPlan{nullptr};             
  std::vector<float*> mVal1DIDCs;            
  std::vector<fftwf_complex*> mCoefficients; 
 
  void calcFourierCoefficientsNaive();
 
  void calcFourierCoefficientsFFTW3();
 
  std::vector<std::vector<float>> inverseFourierTransformNaive() const;
 
  std::vector<std::vector<float>> inverseFourierTransformFFTW3() const;
 
  void normalizeCoefficients()
  {
    std::transform(mFourierCoefficients.mFourierCoefficients.begin(), mFourierCoefficients.mFourierCoefficients.end(), mFourierCoefficients.mFourierCoefficients.begin(), [norm = this->mRangeIDC](auto& val) { return val / norm; });
  };
 
  unsigned int getNMaxCoefficients() const { return this->mRangeIDC / 2 + 1; }
 
  void fftwLoop(const std::vector<float>& idcOneExpanded, const std::vector<unsigned int>& offsetIndex, const unsigned int interval, const unsigned int thread);
 
  ClassDefNV(IDCFourierTransform, 1)
};
 
} // namespace o2::tpc
 
#endif