#include <ProtoDUNEChannelNoiseService.h>

Inheritance diagram for ProtoDUNEChannelNoiseService:

Public Member Functions
	ProtoDUNEChannelNoiseService (fhicl::ParameterSet const &pset)

	ProtoDUNEChannelNoiseService (fhicl::ParameterSet const &pset, art::ActivityRegistry &)

	~ProtoDUNEChannelNoiseService ()

int	addNoise (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, Channel chan, AdcSignalVector &sigs) const

std::ostream &	print (std::ostream &out=std::cout, std::string prefix="") const

void	generateNoise (detinfo::DetectorClocksData const &clockData, float wirelength, float ENOB, float aLowCutoff, AdcSignalVector &noise, TH1 *aNoiseHist) const

Public Member Functions inherited from ChannelNoiseService
virtual	~ChannelNoiseService ()=default

virtual void	newEvent ()

Private Member Functions
void	generateNoise (detinfo::DetectorClocksData const &clockData)

Private Attributes
float	fLowCutoffZ
	low frequency filter cutoff (kHz) for Z (collection) plane More...

float	fLowCutoffU
	low frequency filter cutoff (kHz) for U plane More...

float	fLowCutoffV
	low frequency filter cutoff (kHz) for V plane More...

unsigned int	fNoiseArrayPoints
	number of points in randomly generated noise array More...

bool	fOldNoiseIndex
	Use old selection of noise array index. More...

float	fWhiteNoiseZ
	Level (per freq bin) for white noise for Z. More...

float	fWhiteNoiseU
	Level (per freq bin) for white noise for U. More...

float	fWhiteNoiseV
	Level (per freq bin) for white noise for V. More...

int	fRandomSeed
	Seed for random number service. If absent or zero, use SeedSvc. More...

int	fLogLevel
	Log message level: 0=quiet, 1=init only, 2+=every event. More...

float	fWirelengthZ

float	fWirelengthU

float	fWirelengthV

float	fENOB

AdcSignalVectorVector	fNoiseZ
	noise on each channel for each time for Z (collection) plane More...

AdcSignalVectorVector	fNoiseU
	noise on each channel for each time for U plane More...

AdcSignalVectorVector	fNoiseV
	noise on each channel for each time for V plane More...

TH1 *	fNoiseHistZ
	distribution of noise counts for Z More...

TH1 *	fNoiseHistU
	distribution of noise counts for U More...

TH1 *	fNoiseHistV
	distribution of noise counts for V More...

TH1 *	fNoiseChanHist
	distribution of accessed noise samples More...

CLHEP::HepRandomEngine *	m_pran

Additional Inherited Members
Public Types inherited from ChannelNoiseService
typedef unsigned int	Channel

Detailed Description

Definition at line 25 of file ProtoDUNEChannelNoiseService.h.

Constructor & Destructor Documentation

ProtoDUNEChannelNoiseService::ProtoDUNEChannelNoiseService ( fhicl::ParameterSet const & pset )

Definition at line 30 of file ProtoDUNEChannelNoiseService_service.cc.

 : fRandomSeed(0), fLogLevel(1),
   fNoiseHistZ(nullptr), fNoiseHistU(nullptr), fNoiseHistV(nullptr),
   fNoiseChanHist(nullptr),
   m_pran(nullptr) {
   const string myname = "ProtoDUNEChannelNoiseService::ctor: ";
   fLowCutoffZ        = pset.get<double>("LowCutoffZ");
   fLowCutoffU        = pset.get<double>("LowCutoffU");
   fLowCutoffV        = pset.get<double>("LowCutoffV");
   fWhiteNoiseZ       = pset.get<double>("WhiteNoiseZ");
   fWhiteNoiseU       = pset.get<double>("WhiteNoiseU");
   fWhiteNoiseV       = pset.get<double>("WhiteNoiseV");
   fNoiseArrayPoints  = pset.get<unsigned int>("NoiseArrayPoints");
   fOldNoiseIndex     = pset.get<bool>("OldNoiseIndex");
   fWirelengthZ       = pset.get<double>("WireLengthZ");
   fWirelengthU       = pset.get<double>("WireLengthU");
   fWirelengthV       = pset.get<double>("WireLengthV");
   fENOB              = pset.get<double>("EffectiveNBits");
   bool haveSeed = pset.get_if_present<int>("RandomSeed", fRandomSeed);
   if ( fRandomSeed == 0 ) haveSeed = false;
   pset.get_if_present<int>("LogLevel", fLogLevel);
   fNoiseZ.resize(fNoiseArrayPoints);
   fNoiseU.resize(fNoiseArrayPoints);
   fNoiseV.resize(fNoiseArrayPoints);
   int seed = fRandomSeed;
   art::ServiceHandle<art::TFileService> tfs;
   fNoiseHistZ = tfs->make<TH1F>("znoise", ";Z Noise [ADC counts];", 1000,   -10., 10.);
   fNoiseHistU = tfs->make<TH1F>("unoise", ";U Noise [ADC counts];", 1000,   -10., 10.);
   fNoiseHistV = tfs->make<TH1F>("vnoise", ";V Noise [ADC counts];", 1000,   -10., 10.);
   fNoiseChanHist = tfs->make<TH1F>("NoiseChan", ";Noise channel;", fNoiseArrayPoints, 0, fNoiseArrayPoints);
   // Assign a unique name for the random number engine ProtoDUNEChannelNoiseServiceVIII
   // III = for each instance of this class.
   string rname = "ProtoDUNEChannelNoiseService";
   if ( haveSeed ) {
     if ( fLogLevel > 0 ) cout << myname << "WARNING: Using hardwired seed." << endl;
     m_pran = new HepJamesRandom(seed);
   } else {
     if ( fLogLevel > 0 ) cout << myname << "Using NuRandomService." << endl;
     art::ServiceHandle<NuRandomService> seedSvc;
     m_pran = new HepJamesRandom;
     if ( fLogLevel > 0 ) cout << myname << "    Initial seed: " << m_pran->getSeed() << endl;
     seedSvc->registerEngine(NuRandomService::CLHEPengineSeeder(m_pran), rname);
   }
   if ( fLogLevel > 0 ) cout << myname << "  Registered seed: " << m_pran->getSeed() << endl;
   auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
   for ( unsigned int isam=0; isam<fNoiseArrayPoints; ++isam ) {
     generateNoise(clockData, fWirelengthZ, fENOB, fLowCutoffZ, fNoiseZ[isam], fNoiseHistZ);
     generateNoise(clockData, fWirelengthU, fENOB, fLowCutoffU, fNoiseU[isam], fNoiseHistU);
     generateNoise(clockData, fWirelengthV, fENOB, fLowCutoffV, fNoiseV[isam], fNoiseHistV);
   }
   if ( fLogLevel > 1 ) print() << endl;
 }

ProtoDUNEChannelNoiseService::ProtoDUNEChannelNoiseService	(	fhicl::ParameterSet const &	pset,
		art::ActivityRegistry &
	)

Definition at line 86 of file ProtoDUNEChannelNoiseService_service.cc.

87 : ProtoDUNEChannelNoiseService(pset) { }

ProtoDUNEChannelNoiseService::ProtoDUNEChannelNoiseService

ProtoDUNEChannelNoiseService(fhicl::ParameterSet const &pset)

Definition: ProtoDUNEChannelNoiseService_service.cc:30

ProtoDUNEChannelNoiseService::~ProtoDUNEChannelNoiseService ( )

Definition at line 91 of file ProtoDUNEChannelNoiseService_service.cc.

                                                             {
   const string myname = "ProtoDUNEChannelNoiseService::dtor: ";
   if ( fLogLevel > 0 ) {
     cout << myname << "Deleting random engine with seed " << m_pran->getSeed() << endl;
   }
   delete m_pran;
 }

Member Function Documentation

int ProtoDUNEChannelNoiseService::addNoise	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		Channel	chan,
		AdcSignalVector &	sigs
	)		const

virtual

Implements ChannelNoiseService.

Definition at line 101 of file ProtoDUNEChannelNoiseService_service.cc.

                                                                                       {
   CLHEP::RandFlat flat(*m_pran);
   CLHEP::RandGauss gaus(*m_pran);
   unsigned int noisechan = 0;
   if ( fOldNoiseIndex ) {
     // Keep this strange way of choosing noise channel to be consistent with old results.
     // The relative weights of the first and last channels are 0.5 and 0.6.
     noisechan = nearbyint(flat.fire()*(1.*(fNoiseArrayPoints-1)+0.1));
   } else {
     noisechan = flat.fire()*fNoiseArrayPoints;
     if ( noisechan == fNoiseArrayPoints ) --noisechan;
   }
   fNoiseChanHist->Fill(noisechan);
   art::ServiceHandle<geo::Geometry> geo;
   const geo::View_t view = geo->View(chan);
   for ( unsigned int itck=0; itck<sigs.size(); ++itck ) {
     double tnoise = 0.0;
     double wnoise = 0.0;
     if ( view==geo::kU ) {
       tnoise = fNoiseU[noisechan][itck];
       wnoise = fWhiteNoiseU;
     } else if ( view==geo::kV ) {
       tnoise = fNoiseV[noisechan][itck];
       wnoise = fWhiteNoiseV;
     } else {
       tnoise = fNoiseZ[noisechan][itck];
       wnoise = fWhiteNoiseZ;
     }
     if ( wnoise != 0.0 ) tnoise += wnoise*gaus.fire();
     sigs[itck] += tnoise;
   }
   return 0;
 }

void ProtoDUNEChannelNoiseService::generateNoise	(	detinfo::DetectorClocksData const &	clockData,
		float	wirelength,
		float	ENOB,
		float	aLowCutoff,
		AdcSignalVector &	noise,
		TH1 *	aNoiseHist
	)		const

Definition at line 162 of file ProtoDUNEChannelNoiseService_service.cc.

                                                              {
   const string myname = "ProtoDUNEChannelNoiseService::generateNoise: ";
   if ( fLogLevel > 1 ) {
     cout << myname << "Generating noise." << endl;
     if ( fLogLevel > 2 ) {
       cout << myname << "  Cutoff: " << aLowCutoff << endl;
       cout << myname << "    Seed: " << m_pran->getSeed() << endl;
     }
   }
   // Fetch sampling rate.
   float sampleRate = sampling_rate(clockData);
   // Fetch FFT service and # ticks.
   art::ServiceHandle<util::LArFFT> pfft;
   unsigned int ntick = pfft->FFTSize();
   CLHEP::RandFlat flat(*m_pran);
   // Create noise spectrum in frequency.
   unsigned nbin = ntick/2 + 1;
   std::vector<TComplex> noiseFrequency(nbin, 0.);
   double pval = 0.;
 //  double lofilter = 0.;
   double phase = 0.;
   double rnd[3] = {0.};
   
   ////////////////////////////// MicroBooNE noise model/////////////////////////////////
   // vars
   double params[1] = {0.};
   double fitpar[9] = {0.};
   double wldparams[2] = {0.};
   
   // calculate FFT magnitude of noise from ENOB
   double baseline_noise = std::sqrt(ntick*1.0/12)*std::pow(2, 12 - fENOB);
   // wire length dependence function 
   TF1* _wld_f = new TF1("_wld_f", "[0] + [1]*x", 0.0, 1000);
   // custom poisson  
   TF1* _poisson = new TF1("_poisson", "[0]**(x) * exp(-[0]) / ROOT::Math::tgamma(x+1.)", 0, 30);
   // gain function in MHz
   TF1* _pfn_f1 = new TF1("_pfn_f1", "([0]*1/(x*[8]/2 + 10) + ([1]*exp(-0.5*(((x*[8]/2)-[2])/[3])**2)*exp(-0.5*pow(x*[8]/(2*[4]),[5])))*[6])+[7]", 0.0, (double)ntick);
   
   // set data-driven parameters
   // poisson mean
   params[0] = 3.30762;
   //wire length dependence parameters
   wldparams[0] = 0.395;
   wldparams[1] = 0.001304;
   _wld_f->SetParameters(wldparams);
   double wldValue = _wld_f->Eval(wirelength);
   fitpar[0] = 9.27790e+02;
   fitpar[1] = 1.20284e+07;
   fitpar[2] = 4.93692e+03;
   fitpar[3] = 1.03438e+03;
   fitpar[4] = 2.33306e+02;
   fitpar[5] = 1.36605e+00;
   fitpar[6] = wldValue;
   fitpar[7] = baseline_noise;
   fitpar[8] = ntick;
   _pfn_f1->SetParameters(fitpar);
   _poisson->SetParameters(params);
         
   // width of frequencyBin in kHz
   double binWidth = 1.0/(ntick*sampleRate*1.0e-6);
   for ( unsigned int i=0; i<ntick/2+1; ++i ) {
     //MicroBooNE noise model
     double pfnf1val = _pfn_f1->Eval((double)i*binWidth/1000); //convert to MHz
     // define FFT parameters
     double randomizer = _poisson->GetRandom()/params[0];
     pval = pfnf1val * randomizer;
     // low frequency cutoff     
     //lofilter = 1.0/(1.0+exp(-(i-aLowCutoff/binWidth)/0.5));
     // randomize 10%
     flat.fireArray(2, rnd, 0, 1);
     //pval *= lofilter*(0.9 + 0.2*rnd[0]);
     // random phase angle
     phase = rnd[1]*2.*TMath::Pi();
     TComplex tc(pval*cos(phase),pval*sin(phase));
     noiseFrequency[i] += tc;
   }
   // Obtain time spectrum from frequency spectrum.
   noise.clear();
   noise.resize(ntick,0.0);
   std::vector<double> tmpnoise(noise.size());
   pfft->DoInvFFT(noiseFrequency, tmpnoise);
   noiseFrequency.clear();
   for ( unsigned int itck=0; itck<noise.size(); ++itck ) {
     noise[itck] = /*sqrt(ntick)**/tmpnoise[itck];
   }
   for ( unsigned int itck=0; itck<noise.size(); ++itck ) {
     aNoiseHist->Fill(noise[itck]);
   }
 }

void ProtoDUNEChannelNoiseService::generateNoise ( detinfo::DetectorClocksData const & clockData )

private

Definition at line 256 of file ProtoDUNEChannelNoiseService_service.cc.

                                                                                            {
   fNoiseZ.resize(fNoiseArrayPoints);
   fNoiseU.resize(fNoiseArrayPoints);
   fNoiseV.resize(fNoiseArrayPoints);
   for ( unsigned int inch=0; inch<fNoiseArrayPoints; ++inch ) {
     generateNoise(clockData, fWirelengthZ, fENOB, fLowCutoffZ, fNoiseZ[inch], fNoiseHistZ);
     generateNoise(clockData, fWirelengthU, fENOB, fLowCutoffU, fNoiseZ[inch], fNoiseHistU);
     generateNoise(clockData, fWirelengthV, fENOB, fLowCutoffV, fNoiseZ[inch], fNoiseHistV);
   }
 }

ostream & ProtoDUNEChannelNoiseService::print	(	std::ostream &	out = `std::cout`,
		std::string	prefix = `""`
	)		const

virtual

Implements ChannelNoiseService.

Definition at line 139 of file ProtoDUNEChannelNoiseService_service.cc.

                                                                               {
   out << prefix << "ProtoDUNEChannelNoiseService: " << endl;
   out << prefix << "        LowCutoffZ: " << fLowCutoffZ << endl;
   out << prefix << "        LowCutoffU: " << fLowCutoffU << endl;
   out << prefix << "       WireLengthZ: " << fWirelengthZ  << endl;
   out << prefix << "       WireLengthU: " << fWirelengthU  << endl;
   out << prefix << "       WireLengthV: " << fWirelengthV  << endl;
   out << prefix << "       EffectiveNBits: " << fENOB  << endl;
   out << prefix << "        LowCutoffV: " << fLowCutoffV << endl;
   out << prefix << "  NoiseArrayPoints: " << fNoiseArrayPoints << endl;
   out << prefix << "     OldNoiseIndex: " << fOldNoiseIndex << endl;
   out << prefix << "       WhiteNoiseZ: " << fWhiteNoiseZ << endl;
   out << prefix << "       WhiteNoiseU: " << fWhiteNoiseU << endl;
   out << prefix << "       WhiteNoiseV: " << fWhiteNoiseV << endl;
   out << prefix << "        RandomSeed: " <<  fRandomSeed << endl;
   out << prefix << "          LogLevel: " <<  fLogLevel << endl;
   out << prefix << "  Actual random seed: " << m_pran->getSeed();
   return out;
 }