#include <DigiArapucaSBNDAlg.hh>

Classes
struct	ConfigurationParameters_t

Public Member Functions
	DigiArapucaSBNDAlg (ConfigurationParameters_t const &config)

	~DigiArapucaSBNDAlg ()

double	Baseline ()

void	ConstructWaveform (int ch, sim::SimPhotons const &simphotons, std::vector< short unsigned int > &waveform, std::string pdtype, bool is_daphne, double start_time, unsigned n_samples)

void	ConstructWaveformLite (int ch, sim::SimPhotonsLite const &litesimphotons, std::vector< short unsigned int > &waveform, std::string pdtype, bool is_daphne, double start_time, unsigned n_samples)

Private Member Functions
void	CreatePDWaveform (sim::SimPhotons const &SimPhotons, double t_min, std::vector< double > &wave, std::string pdtype, bool is_daphne)

void	CreatePDWaveformLite (std::map< int, int > const &photonMap, double t_min, std::vector< double > &wave, std::string pdtype, bool is_daphne)

void	SinglePDWaveformCreatorLite (double effT, std::unique_ptr< CLHEP::RandGeneral > &timeHisto, std::vector< double > &wave, std::map< int, int > const &photonMap, double const &t_min, bool is_daphne)

void	SinglePDWaveformCreatorLite (double effT, std::vector< double > &wave, std::map< int, int > const &photonMap, double const &t_min, bool is_daphne)

void	AddSPE (size_t time_bin, std::vector< double > &wave, const std::vector< double > &fWaveformSP, int nphotons)

void	Pulse1PE (std::vector< double > &wave, const double sampling)

void	AddLineNoise (std::vector< double > &wave)

void	AddDarkNoise (std::vector< double > &wave, std::vector< double > &WaveformSP)

double	FindMinimumTime (sim::SimPhotons const &simphotons)

double	FindMinimumTimeLite (std::map< int, int > const &photonMap)

void	CreateSaturation (std::vector< double > &wave)

Private Attributes
const ConfigurationParameters_t	fParams

const double	fSampling

const double	fSampling_Daphne

const double	fXArapucaVUVEff

const double	fXArapucaVISEff

const double	fADCSaturation

CLHEP::HepRandomEngine *	fEngine
	Reference to art-managed random-number engine. More...

CLHEP::RandFlat	fFlatGen

CLHEP::RandPoissonQ	fPoissonQGen

CLHEP::RandGaussQ	fGaussQGen

CLHEP::RandExponential	fExponentialGen

std::unique_ptr < CLHEP::RandGeneral >	fTimeXArapucaVUV

std::unique_ptr < CLHEP::RandGeneral >	fTimeTPB

std::vector< double >	fWaveformSP

std::vector< double >	fWaveformSP_Daphne

std::unordered_map < raw::Channel_t, std::vector < double > >	fFullWaveforms

Detailed Description

Definition at line 41 of file DigiArapucaSBNDAlg.hh.

Constructor & Destructor Documentation

opdet::DigiArapucaSBNDAlg::DigiArapucaSBNDAlg ( ConfigurationParameters_t const & config )

Definition at line 9 of file DigiArapucaSBNDAlg.cc.

     : fParams{config}
     , fSampling(fParams.frequency/ 1000.) //in GHz to cancel with ns
     , fSampling_Daphne(fParams.frequency_Daphne/ 1000.) //in GHz to cancel with ns
     , fXArapucaVUVEff(fParams.XArapucaVUVEff / fParams.larProp->ScintPreScale())
     , fXArapucaVISEff(fParams.XArapucaVISEff / fParams.larProp->ScintPreScale())
     , fADCSaturation(fParams.Baseline + fParams.Saturation * fParams.ADC * fParams.MeanAmplitude)
     , fEngine(fParams.engine)
     , fFlatGen(*fEngine)
     , fPoissonQGen(*fEngine)
     , fGaussQGen(*fEngine)
     , fExponentialGen(*fEngine)
   {
 
     if(fXArapucaVUVEff > 1.0001 || fXArapucaVISEff > 1.0001)
       mf::LogWarning("DigiArapucaSBNDAlg")
         << "Quantum efficiency set in fhicl file "
         << fParams.XArapucaVUVEff << " or " << fParams.XArapucaVISEff
         << " seems to be too large!\n"
         << "Final QE must be equal to or smaller than the scintillation "
         << "pre scale applied at simulation time.\n"
         << "Please check this number (ScintPreScale): "
         << fParams.larProp->ScintPreScale();
 
     std::string fname;
     cet::search_path sp("FW_SEARCH_PATH");
     sp.find_file(fParams.ArapucaDataFile, fname);
     TFile* file = TFile::Open(fname.c_str(), "READ");
     //Note: TPB time now implemented at digitization module for both coated pmts and (x)arapucas
     //OpDetSim/digi_arapuca_sbnd.root updated in sbnd_data (now including the TPB times vector)
 
     // TPB emission time histogram for visible (x)arapucas
     std::vector<double>* timeTPB_p;
     file->GetObject("timeTPB", timeTPB_p);
     fTimeTPB = std::make_unique<CLHEP::RandGeneral>
       (*fEngine, timeTPB_p->data(), timeTPB_p->size());
 
     std::vector<double>* TimeXArapucaVUV_p;
     file->GetObject("TimeXArapucaVUV", TimeXArapucaVUV_p);
     fTimeXArapucaVUV = std::make_unique<CLHEP::RandGeneral>
       (*fEngine, TimeXArapucaVUV_p->data(), TimeXArapucaVUV_p->size());
 
     size_t pulseSize = fParams.PulseLength * fSampling;
     fWaveformSP.resize(pulseSize);
 
     size_t pulseSize_Daphne = fParams.PulseLength * fSampling_Daphne;
     fWaveformSP_Daphne.resize(pulseSize_Daphne);
   
     if(fParams.ArapucaSinglePEmodel) {
       mf::LogDebug("DigiArapucaSBNDAlg") << " using testbench pe response";
       TFile* file =  TFile::Open(fname.c_str(), "READ");
       std::vector<double>* SinglePEVec_40ftCable_Daphne;
       std::vector<double>* SinglePEVec_40ftCable_Apsaia;
       file->GetObject("SinglePEVec_40ftCable_Apsaia", SinglePEVec_40ftCable_Apsaia);
       file->GetObject("SinglePEVec_40ftCable_Daphne", SinglePEVec_40ftCable_Daphne);
       fWaveformSP = *SinglePEVec_40ftCable_Apsaia;
       fWaveformSP_Daphne = *SinglePEVec_40ftCable_Daphne;
     }  
     else{
       mf::LogDebug("DigiArapucaSBNDAlg") << " using ideal pe response";
       Pulse1PE(fWaveformSP,fSampling);
       Pulse1PE(fWaveformSP_Daphne,fSampling_Daphne);
     }
     file->Close();
   } // end constructor

opdet::DigiArapucaSBNDAlg::~DigiArapucaSBNDAlg ( )

Definition at line 75 of file DigiArapucaSBNDAlg.cc.

75 {}

Member Function Documentation

void opdet::DigiArapucaSBNDAlg::AddDarkNoise	(	std::vector< double > &	wave,
		std::vector< double > &	WaveformSP
	)

private

Definition at line 327 of file DigiArapucaSBNDAlg.cc.

   {
     int nCT;
     // Multiply by 10^9 since fDarkNoiseRate is in Hz (conversion from s to ns)
     double mean = 1000000000.0 / fParams.DarkNoiseRate;
     double darkNoiseTime = fExponentialGen.fire(mean);
     while(darkNoiseTime < wave.size()) {
       size_t timeBin = std::round(darkNoiseTime);
       if(fParams.CrossTalk > 0.0 && fFlatGen.fire(1.0) < fParams.CrossTalk) nCT = 2;
       else nCT = 1;
       if(timeBin < wave.size()) AddSPE(timeBin, wave, WaveformSP, nCT);
       // Find next time to add dark noise
       darkNoiseTime += fExponentialGen.fire(mean);
     }//while
   }

void opdet::DigiArapucaSBNDAlg::AddLineNoise ( std::vector< double > & wave )

private

Definition at line 306 of file DigiArapucaSBNDAlg.cc.

   {
     // TODO: after running the profiler I can see that this is where
     // most cycles are being used.  Potentially some improvement could
     // be achieved with the below code but the array dynamic allocation
     // is not helping. The function would need to have it passed.
     // ~icaza
     //
     // double *array = new double[wave.size()]();
     // CLHEP::RandGaussQ::shootArray(fEngine, wave.size(), array, 0, fParams.BaselineRMS);
     // for(size_t i = 0; i<wave.size(); i++) {
     //   wave[i] += array[i];
     // }
     // delete array;
     //
     std::transform(wave.begin(), wave.end(), wave.begin(),
                    [this](double w) -> double {
                      return w + fGaussQGen.fire(0., fParams.BaselineRMS) ; });
   }

void opdet::DigiArapucaSBNDAlg::AddSPE	(	size_t	time_bin,
		std::vector< double > &	wave,
		const std::vector< double > &	fWaveformSP,
		int	nphotons
	)

private

Definition at line 279 of file DigiArapucaSBNDAlg.cc.

   {
     // if(time_bin > wave.size()) return;
     size_t max = time_bin + fWaveformSP.size() < wave.size() ? time_bin + fWaveformSP.size() : wave.size();
     auto min_it = std::next(wave.begin(), time_bin);
     auto max_it = std::next(wave.begin(), max);
 
     double nphotons_aux= nphotons;
     if(fParams.MakeAmpFluctuations) nphotons_aux = fGaussQGen.fire(nphotons, std::sqrt(nphotons) * fParams.AmpFluctuation);
 
     std::transform(min_it, max_it,
                    fWaveformSP.begin(), min_it,
                    [nphotons_aux](double w, double ws) -> double {
                      return w + ws*nphotons_aux  ; });
   }

double opdet::DigiArapucaSBNDAlg::Baseline ( )

inline

Definition at line 77 of file DigiArapucaSBNDAlg.hh.

       {
         return fParams.Baseline;
       }

void opdet::DigiArapucaSBNDAlg::ConstructWaveform	(	int	ch,
		sim::SimPhotons const &	simphotons,
		std::vector< short unsigned int > &	waveform,
		std::string	pdtype,
		bool	is_daphne,
		double	start_time,
		unsigned	n_samples
	)

Definition at line 78 of file DigiArapucaSBNDAlg.cc.

   {
     std::vector<double> waves(n_samples, fParams.Baseline);
     CreatePDWaveform(simphotons, start_time, waves, pdtype,is_daphne);
     waveform = std::vector<short unsigned int> (waves.begin(), waves.end());
   }

void opdet::DigiArapucaSBNDAlg::ConstructWaveformLite	(	int	ch,
		sim::SimPhotonsLite const &	litesimphotons,
		std::vector< short unsigned int > &	waveform,
		std::string	pdtype,
		bool	is_daphne,
		double	start_time,
		unsigned	n_samples
	)

Definition at line 93 of file DigiArapucaSBNDAlg.cc.

   {
     std::vector<double> waves(n_samples, fParams.Baseline);
     std::map<int, int> const& photonMap = litesimphotons.DetectedPhotons;
     CreatePDWaveformLite(photonMap, start_time, waves, pdtype,is_daphne);
     // std::ofstream ofs("True_PE.log",std::ofstream::out | std::ofstream::app);
     // ofs<<ch<<"\t"<<P_truth<<std::endl;
     // ofs.close();
     // P_truth=0;
     waveform = std::vector<short unsigned int> (waves.begin(), waves.end());
   }

void opdet::DigiArapucaSBNDAlg::CreatePDWaveform	(	sim::SimPhotons const &	SimPhotons,
		double	t_min,
		std::vector< double > &	wave,
		std::string	pdtype,
		bool	is_daphne
	)

private

Definition at line 113 of file DigiArapucaSBNDAlg.cc.

   {
     int nCT = 1;
     if(pdtype == "xarapuca_vuv") {
       for(size_t i = 0; i < simphotons.size(); i++) {
         if(fFlatGen.fire(1.0) < fXArapucaVUVEff) {
           double tphoton = (fTimeXArapucaVUV->fire()) + simphotons[i].Time - t_min;
           if(tphoton < 0.) continue; // discard if it didn't made it to the acquisition
           if(fParams.CrossTalk > 0.0 && fFlatGen.fire(1.0) < fParams.CrossTalk) nCT = 2;
           else nCT = 1;
           size_t timeBin = (is_daphne) ? std::floor(tphoton * fSampling_Daphne) : std::floor(tphoton * fSampling);
           if(timeBin < wave.size()) {
             if (!is_daphne) {AddSPE(timeBin, wave, fWaveformSP, nCT);
             }
             else{ AddSPE(timeBin, wave, fWaveformSP_Daphne, nCT);}
           }
         }
       }
     }
     else if(pdtype == "xarapuca_vis") {
       for(size_t i = 0; i < simphotons.size(); i++) {
         if(fFlatGen.fire(1.0) < fXArapucaVISEff) {
           double tphoton = fExponentialGen.fire(fParams.DecayTXArapucaVIS) + simphotons[i].Time - t_min + fTimeTPB->fire();
           if(tphoton < 0.) continue; // discard if it didn't made it to the acquisition
           if(fParams.CrossTalk > 0.0 && fFlatGen.fire(1.0) < fParams.CrossTalk) nCT = 2;
           else nCT = 1;
           size_t timeBin = (is_daphne) ? std::floor(tphoton * fSampling_Daphne) : std::floor(tphoton * fSampling);
           if(timeBin < wave.size()) {
             if (!is_daphne) {AddSPE(timeBin, wave, fWaveformSP, nCT);
             }
             else{ AddSPE(timeBin, wave, fWaveformSP_Daphne, nCT);
             }
           }
         }
       }
     }
     else{
       throw cet::exception("DigiARAPUCASBNDAlg") << "Wrong pdtype: " << pdtype << std::endl;
     }
     if(fParams.BaselineRMS > 0.0) AddLineNoise(wave);
     if(fParams.DarkNoiseRate > 0.0)
     {
       if (!is_daphne) AddDarkNoise(wave,fWaveformSP);
       else            AddDarkNoise(wave,fWaveformSP_Daphne);
     } 
     CreateSaturation(wave);
   }

void opdet::DigiArapucaSBNDAlg::CreatePDWaveformLite	(	std::map< int, int > const &	photonMap,
		double	t_min,
		std::vector< double > &	wave,
		std::string	pdtype,
		bool	is_daphne
	)

private

Definition at line 167 of file DigiArapucaSBNDAlg.cc.

   {
     if(pdtype == "xarapuca_vuv"){
       SinglePDWaveformCreatorLite(fXArapucaVUVEff, fTimeXArapucaVUV, wave, photonMap, t_min,is_daphne);
     }
     else if(pdtype == "xarapuca_vis"){
       // creating the waveforms for xarapuca_vis is different than the rest
       // so there's an overload for that which lacks the timeHisto
       SinglePDWaveformCreatorLite(fXArapucaVISEff, wave, photonMap, t_min,is_daphne);
     }
     else{
       throw cet::exception("DigiARAPUCASBNDAlg") << "Wrong pdtype: " << pdtype << std::endl;
     }
     if(fParams.BaselineRMS > 0.0) AddLineNoise(wave);
     if(fParams.DarkNoiseRate > 0.0)
     {
       if (!is_daphne) AddDarkNoise(wave,fWaveformSP);
       else            AddDarkNoise(wave,fWaveformSP_Daphne);
     } 
 
     CreateSaturation(wave);
   }

void opdet::DigiArapucaSBNDAlg::CreateSaturation ( std::vector< double > & wave )

private

Definition at line 299 of file DigiArapucaSBNDAlg.cc.

   {
     std::replace_if(wave.begin(), wave.end(),
                     [&](auto w){return w > fADCSaturation;}, fADCSaturation);
   }

double opdet::DigiArapucaSBNDAlg::FindMinimumTime ( sim::SimPhotons const & simphotons )

private

Definition at line 344 of file DigiArapucaSBNDAlg.cc.

   {
     double t_min = 1e15;
     for(size_t i = 0; i < simphotons.size(); i++) {
       if(simphotons[i].Time < t_min) t_min = simphotons[i].Time;
     }
     return t_min;
   }

double opdet::DigiArapucaSBNDAlg::FindMinimumTimeLite ( std::map< int, int > const & photonMap )

private

Definition at line 354 of file DigiArapucaSBNDAlg.cc.

   {
     for (auto const& mapMember : photonMap) {
       if(mapMember.second != 0) return (double)mapMember.first;
     }
     return 1e5;
   }

void opdet::DigiArapucaSBNDAlg::Pulse1PE	(	std::vector< double > &	wave,
		const double	sampling
	)

private

Definition at line 269 of file DigiArapucaSBNDAlg.cc.

   {
     double constT1 = fParams.ADC * fParams.MeanAmplitude;
     for(size_t i = 0; i<fWaveformSP.size(); i++) {
       double time = static_cast<double>(i) / sampling;
       if (time < fParams.PeakTime) fWaveformSP[i] = constT1 * std::exp((time - fParams.PeakTime) / fParams.RiseTime);
       else fWaveformSP[i] = constT1 * std::exp(-(time - fParams.PeakTime) / fParams.FallTime);
     }
   }

void opdet::DigiArapucaSBNDAlg::SinglePDWaveformCreatorLite	(	double	effT,
		std::unique_ptr< CLHEP::RandGeneral > &	timeHisto,
		std::vector< double > &	wave,
		std::map< int, int > const &	photonMap,
		double const &	t_min,
		bool	is_daphne
	)

private

Definition at line 196 of file DigiArapucaSBNDAlg.cc.

   {
     // TODO: check that this new approach of not using the last
     // (1-accepted_photons) doesn't introduce some bias
     double meanPhotons;
     size_t acceptedPhotons;
     double tphoton;
     for (auto const& photonMember : photonMap) {
       // TODO: check that this new approach of not using the last
       // (1-accepted_photons) doesn't introduce some bias
       meanPhotons = photonMember.second*effT;
       acceptedPhotons = fPoissonQGen.fire(meanPhotons);
       for(size_t i = 0; i < acceptedPhotons; i++) {
         tphoton = timeHisto->fire();
         tphoton += photonMember.first - t_min;
         if(tphoton < 0.) continue; // discard if it didn't made it to the acquisition
         int nCT=1;
         if(fParams.CrossTalk > 0.0 &&
            fFlatGen.fire(1.0) < fParams.CrossTalk) nCT = 2;
           size_t timeBin = (is_daphne) ? std::floor(tphoton * fSampling_Daphne) : std::floor(tphoton * fSampling);
           if(timeBin < wave.size()) {
             // P_truth=P_truth+nCT;
             if (!is_daphne) {AddSPE(timeBin, wave, fWaveformSP, nCT);
             }
             else{ AddSPE(timeBin, wave, fWaveformSP_Daphne, nCT);}
           }
       }
     }
   }

void opdet::DigiArapucaSBNDAlg::SinglePDWaveformCreatorLite	(	double	effT,
		std::vector< double > &	wave,
		std::map< int, int > const &	photonMap,
		double const &	t_min,
		bool	is_daphne
	)

private

Definition at line 234 of file DigiArapucaSBNDAlg.cc.

   {
     double meanPhotons;
     size_t acceptedPhotons;
     double tphoton;
     int nCT;
     for (auto const& photonMember : photonMap) {
       // TODO: check that this new approach of not using the last
       // (1-accepted_photons) doesn't introduce some bias
       meanPhotons = photonMember.second*effT;
       acceptedPhotons = fPoissonQGen.fire(meanPhotons);
       for(size_t i = 0; i < acceptedPhotons; i++) {
         tphoton = fExponentialGen.fire(fParams.DecayTXArapucaVIS);
         tphoton += photonMember.first - t_min + fTimeTPB->fire();
         if(tphoton < 0.) continue; // discard if it didn't made it to the acquisition
         if(fParams.CrossTalk > 0.0 && fFlatGen.fire(1.0) < fParams.CrossTalk) nCT = 2;
         else nCT = 1;
           size_t timeBin = (is_daphne) ? std::floor(tphoton * fSampling_Daphne) : std::floor(tphoton * fSampling);
           if(timeBin < wave.size()) {
             // P_truth=P_truth+nCT;
             if (!is_daphne) {AddSPE(timeBin, wave, fWaveformSP, nCT);
             }
             else{ AddSPE(timeBin, wave, fWaveformSP_Daphne, nCT);}
           }
       }
     }
   }

Member Data Documentation

const double opdet::DigiArapucaSBNDAlg::fADCSaturation

private

Definition at line 107 of file DigiArapucaSBNDAlg.hh.

CLHEP::HepRandomEngine* opdet::DigiArapucaSBNDAlg::fEngine

private

Reference to art-managed random-number engine.

Definition at line 109 of file DigiArapucaSBNDAlg.hh.

CLHEP::RandExponential opdet::DigiArapucaSBNDAlg::fExponentialGen

private

Definition at line 113 of file DigiArapucaSBNDAlg.hh.

CLHEP::RandFlat opdet::DigiArapucaSBNDAlg::fFlatGen

private

Definition at line 110 of file DigiArapucaSBNDAlg.hh.

std::unordered_map< raw::Channel_t, std::vector<double> > opdet::DigiArapucaSBNDAlg::fFullWaveforms

private

Definition at line 119 of file DigiArapucaSBNDAlg.hh.

CLHEP::RandGaussQ opdet::DigiArapucaSBNDAlg::fGaussQGen

private

Definition at line 112 of file DigiArapucaSBNDAlg.hh.

const ConfigurationParameters_t opdet::DigiArapucaSBNDAlg::fParams

private

Definition at line 101 of file DigiArapucaSBNDAlg.hh.

CLHEP::RandPoissonQ opdet::DigiArapucaSBNDAlg::fPoissonQGen

private

Definition at line 111 of file DigiArapucaSBNDAlg.hh.

const double opdet::DigiArapucaSBNDAlg::fSampling

private

Definition at line 103 of file DigiArapucaSBNDAlg.hh.

const double opdet::DigiArapucaSBNDAlg::fSampling_Daphne

private

Definition at line 104 of file DigiArapucaSBNDAlg.hh.

std::unique_ptr<CLHEP::RandGeneral> opdet::DigiArapucaSBNDAlg::fTimeTPB

private

Definition at line 115 of file DigiArapucaSBNDAlg.hh.

std::unique_ptr<CLHEP::RandGeneral> opdet::DigiArapucaSBNDAlg::fTimeXArapucaVUV

private

Definition at line 114 of file DigiArapucaSBNDAlg.hh.

std::vector<double> opdet::DigiArapucaSBNDAlg::fWaveformSP

private

Definition at line 117 of file DigiArapucaSBNDAlg.hh.

std::vector<double> opdet::DigiArapucaSBNDAlg::fWaveformSP_Daphne

private

Definition at line 118 of file DigiArapucaSBNDAlg.hh.

const double opdet::DigiArapucaSBNDAlg::fXArapucaVISEff

private

Definition at line 106 of file DigiArapucaSBNDAlg.hh.

const double opdet::DigiArapucaSBNDAlg::fXArapucaVUVEff

private

Definition at line 105 of file DigiArapucaSBNDAlg.hh.

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

Classes

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation