#include <DPhaseRealisticNoiseService.h>

Inheritance diagram for DPhaseRealisticNoiseService:

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

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

	~DPhaseRealisticNoiseService ()

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::DetectorPropertiesData const &detProp, std::vector< double > frequencyVector, AdcSignalVector &noise, TH1 *aNoiseHist, double customRandom)

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

virtual void	newEvent ()

Private Member Functions
unsigned int	GetModelSize () const

void	SetModelSize (unsigned int size)

double	GetShift (AdcSignalVector time_vector, int window_length) const

void	Chan2Phase (std::map< Channel, double > &PhaseChannelMap) const

void	generateNoise (detinfo::DetectorPropertiesData const &detProp)

void	mirrorWaveform (AdcSignalVector &noise, int TimeSamples) const

void	importNoiseModel (std::string noiseModel, std::vector< double > &frequencyArrayX, std::vector< double > &frequencyArrayY) const

Private Attributes
std::string	fNoiseModel
	noise model root file More...

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

double	fRandomizeX
	randomization of the average frequency spectrum (on kX or kZ) More...

double	fRandomizeY
	< randomization of the average frequency spectrum (on kY) More...

double	fSmooth

bool	fSetFirst0
	< set the first bin of the frequency array to 0 More...

bool	fSetBaseline
	< Sum baseline model to the data More...

bool	fOldNoiseIndex
	Use old selection of noise array index. 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...

std::vector< double >	fNoiseModelFrequenciesX
	Array storing the frequencies imported from the model in kHz for plane kX (kZ) More...

std::vector< double >	fNoiseModelFrequenciesY
	Array storing the frequencies imported from the model in kHz for plane kY (kZ) More...

AdcSignalVectorVector	fNoiseX
	noise on each channel for each time for X plane More...

AdcSignalVectorVector	fNoiseY
	noise on each channel for each time for Y plane More...

TH1 *	fNoiseHistX
	distribution of noise counts for X More...

TH1 *	fNoiseHistY
	distribution of noise counts for Y More...

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

unsigned int	fModelsize

CLHEP::HepRandomEngine *	m_pran

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

Detailed Description

Definition at line 28 of file DPhaseRealisticNoiseService.h.

Constructor & Destructor Documentation

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

Definition at line 44 of file DPhaseRealisticNoiseService_service.cc.

 : fRandomSeed(0), fLogLevel(1),
   fNoiseHistX(nullptr), fNoiseHistY(nullptr),
   fNoiseChanHist(nullptr),
   m_pran(nullptr) {
   const string myname = "DPhaseRealisticNoiseService::ctor: ";
   fNoiseModel        = pset.get<string>("NoiseModel");
   fRandomizeX        = pset.get<double>("RandomizeX");
   fRandomizeY        = pset.get<double>("RandomizeY");
   fSmooth            = pset.get<double>("Smooth");
   fSetFirst0         = pset.get<bool>("SetFirst0");
   fSetBaseline          = pset.get<bool>("SetBaseline");
   fNoiseArrayPoints  = pset.get<unsigned int>("NoiseArrayPoints");
   fOldNoiseIndex     = pset.get<bool>("OldNoiseIndex");
   bool haveSeed = pset.get_if_present<int>("RandomSeed", fRandomSeed);
   if ( fRandomSeed == 0 ) haveSeed = false;
   pset.get_if_present<int>("LogLevel", fLogLevel);
   fNoiseX.resize(fNoiseArrayPoints);
   fNoiseY.resize(fNoiseArrayPoints);
   int seed = fRandomSeed;
   art::ServiceHandle<art::TFileService> tfs;
   fNoiseHistX = tfs->make<TH1F>("xnoise", ";X Noise [ADC counts];", 1000,   -0.1, 0.1);
   fNoiseHistY = tfs->make<TH1F>("ynoise", ";Y Noise [ADC counts];", 1000,   -0.1, 0.1);
   fNoiseChanHist = tfs->make<TH1F>("NoiseChan", ";Noise channel;", fNoiseArrayPoints, 0, fNoiseArrayPoints);
   // Assign a unique name for the random number engine ExponentialChannelNoiseServiceVIII
   // III = for each instance of this class.
   string rname = "DPhaseRealisticNoiseService";
   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;
 
   importNoiseModel(fNoiseModel, fNoiseModelFrequenciesX, fNoiseModelFrequenciesY);
 
   //pregenerated noise model waveforms based on the realistic noise frequency pattern
   auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
   for ( unsigned int isam=0; isam<fNoiseArrayPoints; ++isam ) {
     generateNoise(detProp, fNoiseModelFrequenciesX,  fNoiseX[isam], fNoiseHistX, fRandomizeX);
     generateNoise(detProp, fNoiseModelFrequenciesY,  fNoiseY[isam], fNoiseHistY, fRandomizeY);
   }
 
   if ( fLogLevel > 1 ) print() << endl;
 } //m_pran(nullptr)

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

Definition at line 98 of file DPhaseRealisticNoiseService_service.cc.

99 : DPhaseRealisticNoiseService(pset) { }

DPhaseRealisticNoiseService::DPhaseRealisticNoiseService

DPhaseRealisticNoiseService(fhicl::ParameterSet const &pset)

Definition: DPhaseRealisticNoiseService_service.cc:44

DPhaseRealisticNoiseService::~DPhaseRealisticNoiseService ( )

Definition at line 103 of file DPhaseRealisticNoiseService_service.cc.

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

Member Function Documentation

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

virtual

Implements ChannelNoiseService.

Definition at line 113 of file DPhaseRealisticNoiseService_service.cc.

                                                                                     {
 
   CLHEP::RandFlat flat(*m_pran);
 
   //define the baseline fluctuations
   //done here because chan is needed to tune the phase
   unsigned int ntick = detProp.NumberTimeSamples();
   double dt = 1./sampling_rate(clockData);
   unsigned int model_tick = GetModelSize();
   std::map<Channel, double>  fPhaseChannelMap;
   Chan2Phase(fPhaseChannelMap);
 
   //define a short osclillating baseline function
   double params[3] = {0.};
   TF1 * _sin = new TF1("_sin", "[0]*TMath::Sin( [1]*x + [2] )", 0.,
                                                              (double)ntick*dt);
   params[0] = 1.; //in adc
   params[1] = (2*TMath::Pi())/(2*model_tick*dt);
   params[2] = fPhaseChannelMap[chan];
 
   _sin->SetParameters(params);
 
   //add one of the pregenerated random noise vectors to the signal vector
   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);
 
   //sigs.resize( detProp.NumberTimeSamples() );
 
   for ( unsigned int itck=0; itck<ntick; ++itck ) {
     double tnoise = 0.0;
     if ( view==geo::kY ) {
       tnoise = fNoiseY[noisechan][itck];
     } else if ( view==geo::kZ ) {
       tnoise = fNoiseX[noisechan][itck];
     } else {
       tnoise = fNoiseX[noisechan][itck];
     }
 
     if(fSetBaseline){
       sigs[itck] += _sin->Eval((double)itck*dt);
     }
 
     sigs[itck] += tnoise;
   }
 
   return 0;
 }

void DPhaseRealisticNoiseService::Chan2Phase ( std::map< Channel, double > & PhaseChannelMap ) const

private

Definition at line 272 of file DPhaseRealisticNoiseService_service.cc.

                                                                                           {
   //create an association between channel and phase.
 
   CLHEP::RandFlat flat(*m_pran);
   art::ServiceHandle<geo::Geometry> geo;
 
   double phase = flat.fire();
   double dph = (TMath::Pi()*0.5)/62; //phase small increment. (May be improved)
 
   for(Channel chan=0; chan< geo->Nchannels() ; chan++){
 
     if(chan % 64 == 0){
       phase = flat.fire();
     }else if(chan % 32){
       dph = -dph;
       phase += dph;
     }else{
       phase += dph;
     }
     PhaseChannelMap[chan] = phase*2*TMath::Pi();
   }//end for
   return;
 }

void DPhaseRealisticNoiseService::generateNoise	(	detinfo::DetectorPropertiesData const &	detProp,
		std::vector< double >	frequencyVector,
		AdcSignalVector &	noise,
		TH1 *	aNoiseHist,
		double	customRandom
	)

Definition at line 375 of file DPhaseRealisticNoiseService_service.cc.

                                                                              {
 
   const string myname = "DPRealisticNoiseService::generateNoise: ";
   if ( fLogLevel > 1 ) {
     cout << myname << "Generating noise." << endl;
   }
 
   CLHEP::RandFlat flat(*m_pran);
   CLHEP::RandGauss gaus(*m_pran);
 
   // rayleigh
   //TF1* _rayleigh = new TF1("_rayleigh", "[0]*( x/([1]*[1]) )*TMath::Exp(-( (x*x)/(2*[1]*[1]) ))", 0, 200);
   //_rayleigh->SetParameter(0, 7.80e+5);
 
   unsigned int inputFreqSize = frequencyVector.size();
   unsigned int pointFFT = 2*(inputFreqSize-1);
 
   art::ServiceHandle<util::LArFFT> pfft;
   pfft->ReinitializeFFT( pointFFT ," ", pfft->FFTFitBins() );
 
   unsigned int model_tick = pfft->FFTSize(); //ticks of the time model
   unsigned nbin = model_tick/2 + 1;
   SetModelSize( model_tick );
 
   std::vector<TComplex> noiseFrequency(nbin, 0.);
   double pval = 0.;
   double phase = 0.;
 
   for ( unsigned int i=0; i<model_tick/2+1; ++i ) {
 
     pval = frequencyVector[i];
 
     //Randomize amplitude (Gaussian)
     pval += customRandom*gaus.fire();
 
     //Randomize phase (Constant)
     phase = flat.fire()*2.*TMath::Pi();
     TComplex tc(pval*cos(phase),pval*sin(phase));
     noiseFrequency[i] += tc;
   }
 
   //Perform the InvFFT of the randomized noise frequency model
   noise.clear();
   noise.resize(model_tick, 0.0);
   std::vector<double> tmpnoise(noise.size());
   pfft->DoInvFFT(noiseFrequency, tmpnoise);
   noiseFrequency.clear();
 
   //model is not normalized
   for ( unsigned int itck=0; itck<noise.size(); ++itck ) {
     noise[itck] = tmpnoise[itck];
   }
 
   //here the signal is mirrored until it match the number of time samples for
   //a given detector geometry
   int ntick = detProp.NumberTimeSamples();
 
   if(noise.size()<(size_t)ntick){
     mirrorWaveform(noise, ntick);
   }
 
   for ( unsigned int itck=0; itck<noise.size(); ++itck ) {
     aNoiseHist->Fill(noise[itck]);
   }
 
   //restore fft with the same number of points as it was before
   pfft->ReinitializeFFT( ntick ," ", pfft->FFTFitBins() );
 
   return;
 }//end generateNoise

void DPhaseRealisticNoiseService::generateNoise ( detinfo::DetectorPropertiesData const & detProp )

private

Definition at line 450 of file DPhaseRealisticNoiseService_service.cc.

                                                                                             {
   fNoiseX.resize(fNoiseArrayPoints);
   fNoiseY.resize(fNoiseArrayPoints);
   for ( unsigned int inch=0; inch<fNoiseArrayPoints; ++inch ) {
     generateNoise(detProp, fNoiseModelFrequenciesX, fNoiseX[inch], fNoiseHistX, fRandomizeX);
     generateNoise(detProp, fNoiseModelFrequenciesY, fNoiseY[inch], fNoiseHistY, fRandomizeY);
   }
 }

unsigned int DPhaseRealisticNoiseService::GetModelSize ( ) const

private

Definition at line 262 of file DPhaseRealisticNoiseService_service.cc.

                                                             {
   return fModelsize;
 }

double DPhaseRealisticNoiseService::GetShift	(	AdcSignalVector	time_vector,
		int	window_length
	)		const

private

Definition at line 298 of file DPhaseRealisticNoiseService_service.cc.

                                                                                {
   //get the shift requeired for mirroring the waveform
 
   int size = time_vector.size();
   double sum=0.;
   for(int i = size - window_length; i< size; i++){ sum+=time_vector[i]; }
   double shift = sum/window_length;
 
   return shift *=2;
 }

void DPhaseRealisticNoiseService::importNoiseModel	(	std::string	noiseModel,
		std::vector< double > &	frequencyArrayX,
		std::vector< double > &	frequencyArrayY
	)		const

private

Definition at line 193 of file DPhaseRealisticNoiseService_service.cc.

                                                                                 {
 
   //import noise model
   art::ServiceHandle<geo::Geometry> geo;
 
   TFile *fin = TFile::Open(noiseModel.c_str(), "READ");
   if(!fin->IsOpen()){
       cout << "ERROR!: Can't open file: " << noiseModel << endl;
       return;
   }
   else{
       cout << "File: " << noiseModel << " successfully opened!" << endl;
   }
   //get the histogram
   TIter next( fin->GetListOfKeys() );
   TKey *key;
   TProfile *inputHist = new TProfile(); //default initialisazion
 
   while( (key = (TKey*)next()) ){
 
     string name( key->GetName() );
     std::string keyName( key->GetClassName() ); //parse char to string
     if( keyName == "TProfile"){
       inputHist = (TProfile*)fin->Get(key->GetName());
     }
     else{
       std::cout << "ERROR! Object: " << keyName << " in file " << noiseModel
                                      << "has not the right format!" << std::endl;
       fin->Close();
       return;
     }
 
     //loop over frequencies: eliminate the last bin of the fft
     //if geometry is 3x1x1dp, then use both models, otherwise just take 0
 
     frequencyArrayX.resize( inputHist->GetNbinsX() );
     frequencyArrayY.resize( inputHist->GetNbinsX() );
 
     for(size_t f=0; f<(size_t)inputHist->GetNbinsX() ; f++){
         if (name.find("_0")!=string::npos){
 
           //NB: use just the 3m strips view for Far Detector sim
           frequencyArrayY.at(f) = inputHist->GetBinContent(f);
           frequencyArrayX.at(f) = inputHist->GetBinContent(f);
         }
         else if(name.find("_1")!=string::npos){
           //frequencyArrayX.at(f) = inputHist->GetBinContent(f);
           continue;
         }
         else{
           cout << "not valid view " << endl;
           continue;
         }
     } //end f loop
   } //end event loop
 
   //one may want to set the first bin of the model to 0
   if(fSetFirst0){
     frequencyArrayX[0]=0;
     frequencyArrayY[0]=0;
   }
 
   fin->Close();
   return;
 }

void DPhaseRealisticNoiseService::mirrorWaveform	(	AdcSignalVector &	noise,
		int	TimeSamples
	)		const

private

Definition at line 312 of file DPhaseRealisticNoiseService_service.cc.

                                                                                {
 
   //extend by mirroring the waveform length from the  model to match
   // the one of the detector
 
   int ArraySize = noise.size();
   int n_window = fSmooth; //tick window to smooth the mirroring
   double shift =0.;
 
   //initialize first shift value from the mean in the time window
   shift = GetShift(noise, n_window);
 
   noise.resize(2*ArraySize);
 
   //do first mirror
   for(int s = ArraySize; s<2*ArraySize; s++){
     if(s>TimeSamples){ break; }
     noise[s] = -noise[2*ArraySize -s - 1 ] + shift;
   }
 
   //if after the first pass the detector time window is not reached, do more...
   if(TimeSamples > 2*ArraySize){
 
     int n = ceil(TimeSamples/ArraySize);
     int n_pass = ceil(log(n)/log(2))+1;
 
     for(int pass=1; pass<n_pass; pass++){
 
       if(noise.size() > (size_t)TimeSamples){ break; } //no need for an extra pass
 
       shift = GetShift(noise, n_window);
 
       noise.resize( pow(2,(pass+1))*ArraySize );
 
       for(int s = pow(2,pass)*ArraySize; s<pow(2,(pass+1))*ArraySize; s++){
         if(s > TimeSamples){ break; }
         noise[s] = -noise[pow(2,(pass+1))*ArraySize -s - 1 ] + shift;
       }
     }//end for pass
   }//end if
 
   //Correct the waveform inclination
   double y, sx = 0, sy = 0, sxy = 0, sx2 = 0, sy2 = 0;
 
   for (size_t s = 0; s < (size_t)TimeSamples; ++s)
   {
       y = noise[s]; sx += s; sy += y; sxy += s*y; sx2 += s*s; sy2 += y*y;
   }
   double ssx = sx2 - ((sx * sx) / TimeSamples);
   double c = sxy - ((sx * sy) / TimeSamples);
   double mx = sx / TimeSamples;
   double my = sy / TimeSamples;
   double b = my - ((c / ssx) * mx);
   double a = c / ssx;
 
   for (size_t s = 0; s < (size_t)TimeSamples; ++s) { noise[s] -= (a*s + b); }
 
   return;
 }//end mirrorWaveform

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

virtual

Implements ChannelNoiseService.

Definition at line 175 of file DPhaseRealisticNoiseService_service.cc.

                                                                              {
 
   out << prefix << "DPhaseRealisticNoiseService:  " <<  endl;
   out << prefix << " Noise model source file:     " << fNoiseModel   << endl;
   out << prefix << " RandomizationX:              " << fRandomizeX   << endl;
   out << prefix << " RandomizationY:              " << fRandomizeY   << endl;
   out << prefix << " Smooth                       " << fSmooth       << endl;
   out << prefix << " fSetFirst0                   " << fSetFirst0    << endl;
   out << prefix << " fSetBaseline                 " << fSetBaseline  << endl;
   out << prefix << " RandomSeed:                  " << fRandomSeed   << endl;
   out << prefix << " LogLevel:                    " << fLogLevel     << endl;
   out << prefix << " Actual random seed:          " << m_pran->getSeed();
 
   return out;
 }

void DPhaseRealisticNoiseService::SetModelSize ( unsigned int size )

private

Definition at line 266 of file DPhaseRealisticNoiseService_service.cc.

                                                                {
   fModelsize = size;
 }

Member Data Documentation

int DPhaseRealisticNoiseService::fLogLevel

private

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

Definition at line 91 of file DPhaseRealisticNoiseService.h.

unsigned int DPhaseRealisticNoiseService::fModelsize

private

Definition at line 109 of file DPhaseRealisticNoiseService.h.

unsigned int DPhaseRealisticNoiseService::fNoiseArrayPoints

private

number of points in randomly generated noise array

Definition at line 83 of file DPhaseRealisticNoiseService.h.

TH1* DPhaseRealisticNoiseService::fNoiseChanHist

private

distribution of accessed noise samples

Definition at line 105 of file DPhaseRealisticNoiseService.h.

TH1* DPhaseRealisticNoiseService::fNoiseHistX

private

distribution of noise counts for X

Definition at line 103 of file DPhaseRealisticNoiseService.h.

TH1* DPhaseRealisticNoiseService::fNoiseHistY

private

distribution of noise counts for Y

Definition at line 104 of file DPhaseRealisticNoiseService.h.

std::string DPhaseRealisticNoiseService::fNoiseModel

private

noise model root file

Definition at line 82 of file DPhaseRealisticNoiseService.h.

std::vector<double> DPhaseRealisticNoiseService::fNoiseModelFrequenciesX

private

Array storing the frequencies imported from the model in kHz for plane kX (kZ)

Definition at line 94 of file DPhaseRealisticNoiseService.h.

std::vector<double> DPhaseRealisticNoiseService::fNoiseModelFrequenciesY

private

Array storing the frequencies imported from the model in kHz for plane kY (kZ)

Definition at line 95 of file DPhaseRealisticNoiseService.h.

AdcSignalVectorVector DPhaseRealisticNoiseService::fNoiseX

private

noise on each channel for each time for X plane

Definition at line 98 of file DPhaseRealisticNoiseService.h.

AdcSignalVectorVector DPhaseRealisticNoiseService::fNoiseY

private

noise on each channel for each time for Y plane

Definition at line 99 of file DPhaseRealisticNoiseService.h.

bool DPhaseRealisticNoiseService::fOldNoiseIndex

private

Use old selection of noise array index.

Definition at line 89 of file DPhaseRealisticNoiseService.h.

double DPhaseRealisticNoiseService::fRandomizeX

private

randomization of the average frequency spectrum (on kX or kZ)

Definition at line 84 of file DPhaseRealisticNoiseService.h.

double DPhaseRealisticNoiseService::fRandomizeY

private

< randomization of the average frequency spectrum (on kY)

Definition at line 85 of file DPhaseRealisticNoiseService.h.

int DPhaseRealisticNoiseService::fRandomSeed

private

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

Definition at line 90 of file DPhaseRealisticNoiseService.h.

bool DPhaseRealisticNoiseService::fSetBaseline

private

< Sum baseline model to the data

Definition at line 88 of file DPhaseRealisticNoiseService.h.

bool DPhaseRealisticNoiseService::fSetFirst0

private

< set the first bin of the frequency array to 0

Definition at line 87 of file DPhaseRealisticNoiseService.h.

double DPhaseRealisticNoiseService::fSmooth

private

Definition at line 86 of file DPhaseRealisticNoiseService.h.

CLHEP::HepRandomEngine* DPhaseRealisticNoiseService::m_pran

private

Definition at line 111 of file DPhaseRealisticNoiseService.h.

The documentation for this class was generated from the following files:

dunesim/dunesim/DetSim/Service/DPhaseRealisticNoiseService.h
dunesim/dunesim/DetSim/Service/DPhaseRealisticNoiseService_service.cc

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation