Inheritance diagram for opdet::OpDetDigitizerDUNEDP:

Public Member Functions
	OpDetDigitizerDUNEDP (fhicl::ParameterSet const &)

void	produce (art::Event &)

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
void	AddPulse (size_t timeBin, int scale, std::vector< double > &waveform, FocusList &fl, double Gain) const

double	getQE (int OpDet)

double	getGain (int OpDet)

void	CheckFHiCLParameters () const

void	CreateSinglePEWaveform ()

void	CreatePDWaveform (sim::SimPhotonsLite const &, opdet::OpDetResponseInterface const &, geo::Geometry const &, std::vector< std::vector< double > > &, std::vector< FocusList > &)

void	CreatePDWaveform (sim::SimPhotons const &, opdet::OpDetResponseInterface const &, geo::Geometry const &, std::vector< std::vector< double > > &, std::vector< FocusList > &)

void	AddLineNoise (std::vector< std::vector< double > > &, const std::vector< FocusList > &fls) const

void	AddDarkNoise (std::vector< std::vector< double > > &, std::vector< FocusList > &fls, int opDet)

unsigned short	CrossTalk () const

std::vector< short >	VectorOfDoublesToVectorOfShorts (std::vector< double > const &) const

std::map< size_t, std::vector< short > >	SplitWaveform (std::vector< short > const &, const FocusList &)

double	GetDriftWindow (detinfo::DetectorPropertiesData const &detProp) const

double	TickToTime (size_t tick) const

size_t	TimeToTick (double time) const

int	PMTSaturationFunction (int)

double	SumOfElements (std::vector< double >)

Private Attributes
std::vector< std::string >	fInputModule

double	fSampleFreq

double	fTimeBegin

double	fTimeEnd

double	fVoltageToADC

double	fLineNoiseRMS

double	fDarkNoiseRate

double	fCrossTalk

double	fPedestal

double	fGain

bool	fDefaultSimWindow

bool	fFullWaveformOutput

size_t	fReadoutWindow

size_t	fPreTrigger

int	fPadding

bool	fDigiTree_SSP_LED

std::vector< double >	t_photon

std::vector< int >	op_photon

TTree *	arvore2

std::unique_ptr< pmtana::AlgoSSPLeadingEdge >	fThreshAlg

std::unique_ptr< CLHEP::RandGauss >	fRandGauss

std::unique_ptr< CLHEP::RandExponential >	fRandExponential

std::unique_ptr< CLHEP::RandFlat >	fRandFlat

double	fQE

bool	fCustomQE

std::vector< double >	fQEVector

bool	fCustomGain

std::vector< double >	fGainVector

double	fGainsError

double	fSPEWidth

std::vector< double >	fSinglePEWaveform

bool	fNegativeSignal

bool	fExportWaveformTree

TTree *	fWaveformTree

TTree *	fRunInfo

int	Run

int	SubRun

int	Event

unsigned int	nSamples

double	SampleSize

std::vector< short * >	adc_value

unsigned int	nOpDet

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 105 of file OpDetDigitizerDUNEDP_module.cc.

Constructor & Destructor Documentation

opdet::OpDetDigitizerDUNEDP::OpDetDigitizerDUNEDP ( fhicl::ParameterSet const & pset )

Definition at line 247 of file OpDetDigitizerDUNEDP_module.cc.

   : EDProducer{pset}
   {
 
     // This module produces (infrastructure piece)
     produces< std::vector< raw::OpDetWaveform > >();
 
     // Read the fcl-file
     fInputModule = pset.get<std::vector<std::string>>("InputModule",{"largeant"});
     fVoltageToADC       = pset.get< double >("VoltageToADC"      );
     fLineNoiseRMS       = pset.get< double >("LineNoiseRMS"      );
     fCrossTalk          = pset.get< double >("CrossTalk"         );
     fPedestal           = pset.get< double  >("Pedestal"          );
     fDefaultSimWindow   = pset.get< bool   >("DefaultSimWindow"  );
     fFullWaveformOutput = pset.get< bool   >("FullWaveformOutput");
     fReadoutWindow      = pset.get< size_t >("ReadoutWindow"     );
     fPreTrigger         = pset.get< size_t >("PreTrigger"        );
     fPadding            = pset.get< int    >("Padding"           );
     fGain               = pset.get< double >("Gain"              );
     fDigiTree_SSP_LED   = pset.get< bool   >("SSP_LED_DigiTree"  );
     fNegativeSignal     = pset.get< bool   >("NegativeSignal"    );
     fCustomQE           = pset.get< bool   >("CustomQEperOpDet",false    );
     fCustomGain         = pset.get< bool   >("CustomGainperOpDet",false    );
     fGainsError         = pset.get< double >("GainsError",0.0    ); //Factor of gain smearing in units of the gain value
     fSPEWidth           = pset.get< double >("SPEWidth",0.0    ); //Smear the SPE amplitude by the SPE width, in units of the gain value.
     fExportWaveformTree = pset.get< bool   >("ExportWaveformTree", false);
 
     if(fCustomQE) fQEVector = pset.get<std::vector<double>>("QEVector");
     if(fCustomGain) fGainVector = pset.get<std::vector<double>>("GainVector");
 
     fThreshAlg = std::make_unique< pmtana::AlgoSSPLeadingEdge >
                    (pset.get< fhicl::ParameterSet >("algo_threshold"));
 
     if(fDigiTree_SSP_LED){
         art::ServiceHandle< art::TFileService > tfs;
         arvore2 = tfs->make<TTree>("PhotonData", "Photon_analysis");
         arvore2->Branch("photon_opCh",&op_photon);
         arvore2->Branch("photon_pulse",&t_photon);    
     }
 //
     art::ServiceHandle<OpDigiProperties> odp;
 
     // Obtaining parameters from the DetectorClocksService
     auto const clockData = art::ServiceHandle< detinfo::DetectorClocksService const>()->DataForJob();
     fSampleFreq    = clockData.OpticalClock().Frequency();
     fDarkNoiseRate =  odp->DarkRate();
     fQE = odp->QE();
 
 //    fSampleFreq = odp->SampleFreq(); //Sample Freq in MHz
 
     if (fDefaultSimWindow)
     {
       auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);
       // Assume the readout starts at -1 drift window
       fTimeBegin = -1*GetDriftWindow(detProp);
 
       // Take the TPC readout window size and convert
       // to us with the electronics clock frequency
       fTimeEnd   = detProp.ReadOutWindowSize() / clockData.TPCClock().Frequency();
 
       fPreTrigger = 0;//Since we are using negative times as the pretrigger, PreTrigger is zero.
       //fPreTrigger=GetDriftWindow()*clockData.OpticalClock().Frequency();
       fReadoutWindow = (fTimeEnd- fTimeBegin)*clockData.OpticalClock().Frequency();
     }
     else
     {
       fTimeBegin = odp->TimeBegin();
       fTimeEnd   = odp->TimeEnd();
 
     }
 
     CheckFHiCLParameters();
 
     // Initializing random number engines
     unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
     auto& engine = createEngine(seed);
     fRandGauss       = std::make_unique< CLHEP::RandGauss       >(engine);
     fRandExponential = std::make_unique< CLHEP::RandExponential >(engine);
     fRandFlat        = std::make_unique< CLHEP::RandFlat        >(engine);
 
 
     art::ServiceHandle< geo::Geometry > geometry;
     nOpDet=geometry->NOpDets();
 
     if(fGainsError!=0)
     {
       std::cout << "Smearing gains vector by " << fGainsError << " percent. " << std::endl;
       std::vector<double> newGains;
       for (unsigned int pm=0; pm<nOpDet; pm++) newGains.push_back(fRandGauss->fire(getGain(pm), fGainsError*getGain(pm)));
       std::cout << "Old gains vector: "; for (unsigned int pm=0; pm<nOpDet; pm++) std::cout << getGain(pm) << " "; std::cout << std::endl;
       fCustomGain=true; fGainVector=newGains;
       std::cout << "New gains vector: "; for (unsigned int pm=0; pm<nOpDet; pm++) std::cout << getGain(pm) << " "; std::cout << std::endl;
     }
 
     fSinglePEWaveform = odp->SinglePEWaveform();
     if(fSampleFreq!=250)
     {
        int rebin=(int)(1.0*250/fSampleFreq);
        std::vector <double> aux;
 
        std::cout << "\tbefore: SinglePEWaveform: " << fSinglePEWaveform.size() <<"bins, rebin: "<< rebin<<std::endl;
         if (rebin<(int)fSinglePEWaveform.size())
         {
            for(unsigned int i=0; i<fSinglePEWaveform.size();i++)
            {
               if(i%rebin==0) aux.push_back(fSinglePEWaveform[i]);
               else aux[aux.size()-1]+=fSinglePEWaveform[i];
            }
         }
         else
         {
            aux.resize(1);aux[0]=1.0;
         }
         fSinglePEWaveform=aux;
     }
     double normalization =0;
     for(unsigned int i=0; i<fSinglePEWaveform.size();i++) normalization+=fSinglePEWaveform[i];
     std::cout << "Generating waveforms of " << fTimeEnd - fTimeBegin << "us = "<< (fTimeEnd - fTimeBegin)*fSampleFreq <<" Samples"<< std::endl;
     std::cout << "\tTimeBegin" << fTimeBegin <<" "<< std::endl;
     std::cout << "\tfTimeEnd" << fTimeEnd <<" "<< std::endl;
     std::cout << "\tSampleFreq" << fSampleFreq <<" "<< std::endl;
     std::cout << "\tReadoutWindow" << fReadoutWindow <<" "<< std::endl;
     std::cout << "\tfPreTrigger" << fPreTrigger <<" "<< std::endl;
     std::cout << "\fSinglePEWaveform: " << fSinglePEWaveform.size() <<"bins, norm: "<< normalization<<std::endl;
 
     nSamples = (fTimeEnd - fTimeBegin)*fSampleFreq;
     SampleSize = 1000.0/fSampleFreq;
 
     art::ServiceHandle<art::TFileService> tfs;
 
     if(fExportWaveformTree)
     {
       adc_value.resize(nOpDet);
       for (unsigned int i=0; i<nOpDet; i++)adc_value[i] = (short*)malloc(sizeof(short)*nSamples);
 
       double *vaux; vaux = (double*)malloc(sizeof(double)*nOpDet);
       fRunInfo = tfs->make<TTree>("RunInfo","MonteCarlo Run Info");
       fRunInfo->Branch("SampleSize"   , &SampleSize   , "SampleSize/D"     );
       fRunInfo->Branch("nOpDet"        , &nOpDet   , Form("nOpDet/I")    );
       fRunInfo->Branch("Gains"        , vaux   , Form("Gains[nOpDet]/D")    );
       for (unsigned int i=0; i<nOpDet; i++) vaux[i]=getGain(i);
       for (unsigned int i=0; i<nOpDet; i++) std::cout << " Gain [ " << i << " ] = " << vaux[i] << " ADCxticks/PE" << std::endl;
       fRunInfo->Fill();
 
 
       fWaveformTree = tfs->make<TTree>("WaveformTree","Waveforms Tree");
       fWaveformTree->Branch("Run"       , &Run       , "Run/I"       );
       fWaveformTree->Branch("SubRun"    , &SubRun    , "SubRun/I"    );
       fWaveformTree->Branch("Event"     , &Event     , "Event/I"     );
       fWaveformTree->Branch("NSamples"     , &nSamples     , "NSamples/I"     );
       for (unsigned int i=0; i<nOpDet; i++) fWaveformTree->Branch(Form("adc_channel_%i",i), adc_value[i] , Form("adc_value_%i[NSamples]/S",i)     );
 
     }
 
   }

Member Function Documentation

void opdet::OpDetDigitizerDUNEDP::AddDarkNoise	(	std::vector< std::vector< double > > &	waveforms,
		std::vector< FocusList > &	fls,
		int	opDet
	)

private

Definition at line 797 of file OpDetDigitizerDUNEDP_module.cc.

   {
     int i = 0;
     for (auto& waveform : waveforms)
     {
       // Multiply by 10^6 since fDarkNoiseRate is in Hz
       double darkNoiseTime = static_cast< double >(fRandExponential->
                             fire(1.0/fDarkNoiseRate)*1000000.0) + fTimeBegin;
       while (darkNoiseTime < fTimeEnd)
       {
         size_t timeBin = TimeToTick(darkNoiseTime);
         AddPulse(timeBin, CrossTalk(), waveform, fls[i], getGain(optDet) );
         // Find next time to simulate a single PE pulse
         darkNoiseTime += static_cast< double >
                         (fRandExponential->fire(1.0/fDarkNoiseRate)*1000000.0);
       }
 
       ++i;
     }
   }

void opdet::OpDetDigitizerDUNEDP::AddLineNoise	(	std::vector< std::vector< double > > &	waveforms,
		const std::vector< FocusList > &	fls
	)		const

private

Definition at line 780 of file OpDetDigitizerDUNEDP_module.cc.

   {
     int i = 0;
     for(auto& waveform : waveforms){
       for(unsigned int j = 0; j < fls[i].ranges.size(); ++j){
         const std::pair<int, int>& p = fls[i].ranges[j];
         for(int k = p.first; k <= p.second; ++k){
           waveform[k] += fRandGauss->fire(0, fLineNoiseRMS);
         }
       }
 
       ++i;
     }
   }

void opdet::OpDetDigitizerDUNEDP::AddPulse	(	size_t	timeBin,
		int	scale,
		std::vector< double > &	waveform,
		FocusList &	fl,
		double	Gain
	)		const

private

Definition at line 641 of file OpDetDigitizerDUNEDP_module.cc.

   {
         if(timeBin>waveform.size()) return; //photon out of time
 
     size_t pulseLength = fSinglePEWaveform.size()+1;
     if ((timeBin + fSinglePEWaveform.size()) > waveform.size())
       pulseLength = (waveform.size() - timeBin);
 
     fl.AddRange(timeBin, timeBin+pulseLength-1);
 
 
     double interbinallignment = fRandFlat->fire(1.0);//shifting the deposited charge among consecutive bins.
     double customGain = fRandGauss->fire(Gain, fSPEWidth*Gain); //smear the SPE amplitude by the SPEwidth
     // Adding a pulse to the waveform
     for (size_t tick = 0; tick != pulseLength; ++tick)
     {
       double bincharge;
       if (tick==0) bincharge=interbinallignment*scale*customGain*fSinglePEWaveform[tick];
       else bincharge = interbinallignment*scale*customGain*fSinglePEWaveform[tick] + (1-interbinallignment)*scale*customGain*fSinglePEWaveform[tick-1];
 
       if(!fNegativeSignal) waveform[timeBin + tick] += bincharge;
       else waveform[timeBin + tick] -= (double)bincharge;
 
     }
   }

void opdet::OpDetDigitizerDUNEDP::CheckFHiCLParameters ( ) const

private

Definition at line 932 of file OpDetDigitizerDUNEDP_module.cc.

   {
 
     // Are all these logic errors?
 
     if (fLineNoiseRMS < 0.0)
       throw art::Exception(art::errors::LogicError)
                                  << "fLineNoiseRMS: " << fLineNoiseRMS << '\n'
                                  << "Line noise RMS should be non-negative!\n";
 
     if (fDarkNoiseRate < 0.0)
       throw art::Exception(art::errors::LogicError)
                                 << "fDarkNoiseRate: " << fDarkNoiseRate << '\n'
                                 << "Dark noise rate should be non-negative!\n";
 
     if (fPreTrigger >= fReadoutWindow)
       throw art::Exception(art::errors::LogicError)
                << "PreTrigger: "    << fPreTrigger    << " and "
                << "ReadoutWindow: " << fReadoutWindow << '\n'
                << "Pretrigger window has to be shorter than readout window!\n";
 
     if (fTimeBegin >= fTimeEnd)
       throw art::Exception(art::errors::LogicError)
                                  << "TimeBegin: " << fTimeBegin << " and "
                                  << "TimeEnd: "   << fTimeEnd   << '\n'
                                  << "TimeBegin should be less than TimeEnd!\n";
 
   }

void opdet::OpDetDigitizerDUNEDP::CreatePDWaveform	(	sim::SimPhotonsLite const &	litePhotons,
		opdet::OpDetResponseInterface const &	odResponse,
		geo::Geometry const &	geometry,
		std::vector< std::vector< double > > &	pdWaveforms,
		std::vector< FocusList > &	fls
	)

private

Definition at line 671 of file OpDetDigitizerDUNEDP_module.cc.

   {
 
     unsigned int const opDet = litePhotons.OpChannel;
     // This is int because otherwise detectedLite doesn't work
     int readoutChannel;
     // For a group of photons arriving at the same time this is a map
     // of < arrival time (in ns), number of photons >
     std::map< int, int > const& photonsMap = litePhotons.DetectedPhotons;
 
     int counter=0;
     // For every pair of (arrival time, number of photons) in the map:
 
     for (auto const& pulse : photonsMap)
     {
       // Converting ns to us
       double photonTime = static_cast< double >(pulse.first)/1000.0;
 
       int NumberOfPEs = PMTSaturationFunction(pulse.second);
 //      std::cout << "Adding " <<NumberOfPEs << "PEs @ " << photonTime <<  " in ch" << opDet << std::endl;
       for (int i = 0; i < NumberOfPEs; ++i)
       {
         if ((photonTime >= fTimeBegin) && (photonTime < fTimeEnd))
         {
           // Sample a random subset according to QE
           if (CLHEP::RandFlat::shoot(1.0) <getQE(opDet))
           {
             odResponse.detectedLite(opDet, readoutChannel);
             unsigned int hardwareChannel = geometry.HardwareChannelFromOpChannel(readoutChannel);
             // Convert the time of the pulse to ticks
             size_t timeBin = TimeToTick(photonTime);
             // Add 1 pulse to the waveform
             AddPulse(timeBin, CrossTalk(), pdWaveforms.at(hardwareChannel), fls[hardwareChannel], getGain(opDet));
             counter++;
             unsigned int opChannel = geometry.OpChannel(opDet, hardwareChannel);
             if(fDigiTree_SSP_LED){
                 op_photon.emplace_back(opChannel);
                 t_photon.emplace_back(photonTime);
             }
           }
           //else std::cout << "photon not detected "<< std::endl;
         }
       }
     }
 //   std::cout << "Created waveform for channel " << opDet << " with " << counter << " photons." << std::endl;
   
   }

void opdet::OpDetDigitizerDUNEDP::CreatePDWaveform	(	sim::SimPhotons const &	Photons,
		opdet::OpDetResponseInterface const &	odResponse,
		geo::Geometry const &	geometry,
		std::vector< std::vector< double > > &	pdWaveforms,
		std::vector< FocusList > &	fls
	)

private

Definition at line 725 of file OpDetDigitizerDUNEDP_module.cc.

   {
 
     unsigned int const opDet = Photons.OpChannel();
     int readoutChannel;
     // For a group of photons arriving at the same time this is a map
     // of < arrival time (in ns), number of photons >
     std::map< int, int > photonsMap;
     for(unsigned int i=0; i<Photons.size();i++)photonsMap[Photons[i].Time]=0;
     for(unsigned int i=0; i<Photons.size();i++)photonsMap[Photons[i].Time]++;
     int counter=0;
     // For every pair of (arrival time, number of photons) in the map:
 
     for (auto const& pulse : photonsMap)
     {
       // Converting ns to us
       double photonTime = static_cast< double >(pulse.first)/1000.0;
 
       int NumberOfPEs = PMTSaturationFunction(pulse.second);
 //      std::cout << "Adding " <<NumberOfPEs << "PEs @ " << photonTime <<  " in ch" << opDet << std::endl;
       for (int i = 0; i < NumberOfPEs; ++i)
       {
         if ((photonTime >= fTimeBegin) && (photonTime < fTimeEnd))
         {
           // Sample a random subset according to QE
           if (CLHEP::RandFlat::shoot(1.0) <getQE(opDet))
           {
             odResponse.detectedLite(opDet, readoutChannel);
             unsigned int hardwareChannel = geometry.HardwareChannelFromOpChannel(readoutChannel);
             // Convert the time of the pulse to ticks
             size_t timeBin = TimeToTick(photonTime);
             // Add 1 pulse to the waveform
             AddPulse(timeBin, CrossTalk(), pdWaveforms.at(hardwareChannel), fls[hardwareChannel], getGain(opDet));
             counter++;
             unsigned int opChannel = geometry.OpChannel(opDet, hardwareChannel);
             if(fDigiTree_SSP_LED){
                 op_photon.emplace_back(opChannel);
                 t_photon.emplace_back(photonTime);
             }
           }
           //else std::cout << "photon not detected "<< std::endl;
         }
       }
     }
 //   std::cout << "Created waveform for channel " << opDet << " with " << counter << " photons." << std::endl;
   
   }

void opdet::OpDetDigitizerDUNEDP::CreateSinglePEWaveform ( )

private

unsigned short opdet::OpDetDigitizerDUNEDP::CrossTalk ( ) const

private

Definition at line 820 of file OpDetDigitizerDUNEDP_module.cc.

   {
     // Sometimes this should produce 3 or more PEs (not implemented)
     if      (fCrossTalk <= 0.0)                 return 1;
     else if (fRandFlat->fire(1.0) > fCrossTalk) return 1;
     else                                        return 2;
   }

double opdet::OpDetDigitizerDUNEDP::GetDriftWindow ( detinfo::DetectorPropertiesData const & detProp ) const

private

Definition at line 890 of file OpDetDigitizerDUNEDP_module.cc.

   {
 
     double driftWindow;
 
     double maxDrift = 0.0;
     for (geo::TPCGeo const& tpc :
            art::ServiceHandle< geo::Geometry >()->IterateTPCs())
       if (maxDrift < tpc.DriftDistance()) maxDrift = tpc.DriftDistance();
 
     driftWindow = maxDrift / detProp.DriftVelocity();
 
     return driftWindow;
 
   }

double opdet::OpDetDigitizerDUNEDP::getGain ( int OpDet )

private

Definition at line 974 of file OpDetDigitizerDUNEDP_module.cc.

   {
     if(!fCustomGain) return fGain;
     else return fGainVector[OpDet];
   }

double opdet::OpDetDigitizerDUNEDP::getQE ( int OpDet )

private

Definition at line 968 of file OpDetDigitizerDUNEDP_module.cc.

   {
     if(!fCustomQE) return fQE;
     else return fQEVector[OpDet];
   }

int opdet::OpDetDigitizerDUNEDP::PMTSaturationFunction ( int photons )

private

Definition at line 961 of file OpDetDigitizerDUNEDP_module.cc.

961 { return photons;}

void opdet::OpDetDigitizerDUNEDP::produce ( art::Event & evt )

virtual

Implements art::EDProducer.

Definition at line 404 of file OpDetDigitizerDUNEDP_module.cc.

   {
     if(fExportWaveformTree)
     {
       Run    = evt.run();
       SubRun = evt.subRun();
       Event  = evt.event();
     }
     if(fDigiTree_SSP_LED) art::ServiceHandle< art::TFileService > tfs;
 
     // A pointer that will store produced OpDetWaveforms
     std::unique_ptr< std::vector< raw::OpDetWaveform > >
     pulseVecPtr(std::make_unique< std::vector< raw::OpDetWaveform > >());
 
     art::ServiceHandle< sim::LArG4Parameters > lgp;
     bool fUseLitePhotons = lgp->UseLitePhotons();
 
     if (!fUseLitePhotons)
       throw art::Exception(art::errors::UnimplementedFeature)
         << "Sorry, but for now only Lite Photon digitization is implemented!"
         << '\n';
 
     // Geometry service
     art::ServiceHandle< geo::Geometry > geometry;
 
     // Service for determining optical detector responses
     art::ServiceHandle< opdet::OpDetResponseInterface > odResponse;
 
 //    std::vector<art::Handle< std::vector< sim::SimPhotonsLite > > > litePhotonHandle;
 
     int modulecounter=0;
 
 
 //     unsigned int nChannelsPerOpDet=1;
 
      unsigned int nChannelsPerOpDet = geometry->NOpHardwareChannels(nOpDet);
 
 
    // This vector stores waveforms created for each optical channel
     
    for(unsigned int opDet=0; opDet<nOpDet; opDet++)
    {
 
     std::vector< std::vector< double > > pdWaveforms(nChannelsPerOpDet, std::vector< double >(nSamples, static_cast< double >(fPedestal)));
     std::vector<std::vector<FocusList>> fls(nOpDet,std::vector<FocusList>(nChannelsPerOpDet, FocusList(nSamples, fPadding)));
 
     for(auto mod : fInputModule){
       // For every module collection
       //std::cout << "ANALYZING " << mod<<std::endl;                                       
       modulecounter++;
 
       art::Handle< std::vector< sim::SimPhotonsLite > > litePhotonHandle;
       if(fUseLitePhotons) litePhotonHandle = evt.getHandle< std::vector< sim::SimPhotonsLite > >(mod);
 
       art::Handle< std::vector< sim::SimPhotons > > PhotonHandle;
       if(!fUseLitePhotons) PhotonHandle = evt.getHandle< std::vector< sim::SimPhotons > >(mod);
 
       if(fUseLitePhotons)
       {
 
        // For every optical detector:
        for (auto const& litePhotons : (*litePhotonHandle))
        {
           // OpChannel in SimPhotonsLite is actually the photon detector number
           //unsigned int opDet = litePhotons.OpChannel;
          if(opDet == (unsigned)litePhotons.OpChannel)
          {
           // Get number of channels in this optical detector
           //unsigned int nChannelsPerOpDet = geometry->NOpHardwareChannels(opDet);
           // it is one by default in the dual phase geometry, but let's keep it to have compatible functions with single phase.
 
           CreatePDWaveform(litePhotons, *odResponse, *geometry, pdWaveforms, fls[opDet]);
           if((unsigned)modulecounter<fInputModule.size()) continue;//==fInputModule.size()
 
           if(fExportWaveformTree) for (unsigned int hardwareChannel = 0;
              hardwareChannel < nChannelsPerOpDet; ++hardwareChannel) fls[opDet][hardwareChannel].AddRange(0,SampleSize);
 
 
           // Generate dark noise
           if (fDarkNoiseRate > 0.0) AddDarkNoise(pdWaveforms, fls[opDet], opDet);
 
           // Uncomment to undo the effect of FocusLists. Replaces the accumulated
           // lists with ones asserting we need to look at the whole trace.
           // for(FocusList& fl: fls){
           //        fl.ranges.clear();
           //        fl.ranges.emplace_back(0, nSamples-1);
           // }
 
           // Vary the pedestal
 
           if (fLineNoiseRMS > 0.0)  AddLineNoise(pdWaveforms, fls[opDet]);
 
           // Loop over all the created waveforms, split them into shorter
           // waveforms and use them to initialize OpDetWaveforms
 
           for (unsigned int hardwareChannel = 0;
              hardwareChannel < nChannelsPerOpDet; ++hardwareChannel)
           {
             for(const std::pair<int, int>& p: fls[opDet][hardwareChannel].ranges){
             // It's a shame we copy here. We could actually avoid by making the
             // functions below take a begin()/end() pair.
               const std::vector<double> sub(pdWaveforms[hardwareChannel].begin()+p.first,
                                         pdWaveforms[hardwareChannel].begin()+p.second+1);
 
               std::vector< short > waveformOfShorts =
                  VectorOfDoublesToVectorOfShorts(sub);
               //std::cout << "waveformOfShorts " << waveformOfShorts.size()<< std::endl;
 
               std::map< size_t, std::vector < short > > mapTickWaveform =
                 (!fFullWaveformOutput) ?
                 SplitWaveform(waveformOfShorts, fls[opDet][hardwareChannel]) :
                 std::map< size_t, std::vector< short > >{ std::make_pair(0,
                                                                        waveformOfShorts) };
 
               //std::cout << "mapTickWaveform " << mapTickWaveform.size()<< std::endl;
               unsigned int opChannel = geometry->OpChannel(opDet, hardwareChannel);
               for (auto const& pairTickWaveform : mapTickWaveform)
               {
                 double timeStamp =
                 static_cast< double >(TickToTime(pairTickWaveform.first+p.first));
                 //std::cout << "\tp " << p<< std::endl;
                 //std::cout << "\ttimeStamp " << timeStamp<< ", pairTickWaveform.first " << pairTickWaveform.first <<", p.first " << p.first<< std::endl;
                 // std::cout << "\tpairTickWaveform.second.size() " << pairTickWaveform.second.size()<< std::endl;
 
                 raw::OpDetWaveform adcVec(timeStamp, opChannel,
                                           pairTickWaveform.second.size());
                 int counter=0;
                 for (short const& value : pairTickWaveform.second)
                 {
                   //std::cout <<"\t\tvalue " << value<< std::endl;
                   adcVec.emplace_back(value); counter++;
                 } 
                 //std::cout <<"\t\tcounter " << counter<< std::endl;
            
                 pulseVecPtr->emplace_back(std::move(adcVec));
 
               }//endloop per 
             }//endloop per Focus List <fls>
           }//endloop per Hardware channel
          }//endif pmt 
         }//endloop per OpDet
       }//endif litephotons
       else //if SimPhotons
       {  
        // For every optical detector:
        for (auto const& Photons : (*PhotonHandle))
        {
          if(opDet == (unsigned)Photons.OpChannel())
          {
           CreatePDWaveform(Photons, *odResponse, *geometry, pdWaveforms, fls[opDet]);
           if((unsigned)modulecounter<fInputModule.size()) continue;//==fInputModule.size()
 
           if (fDarkNoiseRate > 0.0) AddDarkNoise(pdWaveforms, fls[opDet], opDet);
 
           if (fLineNoiseRMS > 0.0)  AddLineNoise(pdWaveforms, fls[opDet]);
 
           // Loop over all the created waveforms, split them into shorter
           // waveforms and use them to initialize OpDetWaveforms
 
           for (unsigned int hardwareChannel = 0;
              hardwareChannel < nChannelsPerOpDet; ++hardwareChannel)
           {
             for(const std::pair<int, int>& p: fls[opDet][hardwareChannel].ranges){
             // It's a shame we copy here. We could actually avoid by making the
             // functions below take a begin()/end() pair.
               const std::vector<double> sub(pdWaveforms[hardwareChannel].begin()+p.first,
                                         pdWaveforms[hardwareChannel].begin()+p.second+1);
 
               std::vector< short > waveformOfShorts =
                  VectorOfDoublesToVectorOfShorts(sub);
 
               std::map< size_t, std::vector < short > > mapTickWaveform =
                 (!fFullWaveformOutput) ?
                 SplitWaveform(waveformOfShorts, fls[opDet][hardwareChannel]) :
                 std::map< size_t, std::vector< short > >{ std::make_pair(0,
                                                                        waveformOfShorts) };
 
               //std::cout << "mapTickWaveform " << mapTickWaveform.size()<< std::endl;
               unsigned int opChannel = geometry->OpChannel(opDet, hardwareChannel);
               for (auto const& pairTickWaveform : mapTickWaveform)
               {
                 double timeStamp =
                 static_cast< double >(TickToTime(pairTickWaveform.first+p.first));
                 //std::cout << "\tp " << p<< std::endl;
                 //std::cout << "\ttimeStamp " << timeStamp<< ", pairTickWaveform.first " << pairTickWaveform.first <<", p.first " << p.first<< std::endl;
                 // std::cout << "\tpairTickWaveform.second.size() " << pairTickWaveform.second.size()<< std::endl;
 
                 raw::OpDetWaveform adcVec(timeStamp, opChannel,
                                           pairTickWaveform.second.size());
                 int counter=0;
                 for (short const& value : pairTickWaveform.second)
                 {
                   //std::cout <<"\t\tvalue " << value<< std::endl;
                   adcVec.emplace_back(value); counter++;
                 } 
                 //std::cout <<"\t\tcounter " << counter<< std::endl;
            
                 pulseVecPtr->emplace_back(std::move(adcVec));
               }//endloop per 
             }//endloop per Focus List <fls>
           }//endloop per Hardware channel
          }//endif pmt 
         }//endloop per OpDet
       } //endif litePhotons
       if(fUseLitePhotons) litePhotonHandle.clear();
       if(!fUseLitePhotons) PhotonHandle.clear();
     }//endloop per Input Module (S1 and S2 light)
 
    }//endloop per pmt
     if(fDigiTree_SSP_LED)
     {
       arvore2->Fill();
       t_photon.clear();
       op_photon.clear();
     }
 
     if(fExportWaveformTree)
     {
        for(unsigned int j=0; j < pulseVecPtr->size(); j++)
        {
          unsigned int opDet= pulseVecPtr->at(j).ChannelNumber();
 //         std::cout <<" one - pmt - " <<  j << "\t" << opDet  << "\t"<< pulseVecPtr->at(j).Waveform().size() << " " << pulseVecPtr->at(j).Waveform()[0] <<  std::endl;
          for (unsigned int i = 0; i< pulseVecPtr->at(j).Waveform().size() ; i++) adc_value[opDet][i] = pulseVecPtr->at(j).Waveform()[i];
        }
        fWaveformTree->Fill();
 
     }
 
     // Push the OpDetWaveforms into the event
     evt.put(std::move(pulseVecPtr));
 
 
 
 
   }

std::map< size_t, std::vector< short > > opdet::OpDetDigitizerDUNEDP::SplitWaveform	(	std::vector< short > const &	waveform,
		const FocusList &	fls
	)

private

Definition at line 849 of file OpDetDigitizerDUNEDP_module.cc.

   {
 
     std::map< size_t, std::vector< short > > mapTickWaveform;
 
     ::pmtana::PedestalMean_t  ped_mean (waveform.size(),0);
     ::pmtana::PedestalSigma_t ped_sigma(waveform.size(),0);
 
 
     fThreshAlg->Reconstruct(waveform,ped_mean,ped_sigma);
 
     std::vector< pmtana::pulse_param > pulses;
     for (size_t pulseCounter = 0; pulseCounter < fThreshAlg->GetNPulse();
                                                           ++pulseCounter)
       pulses.emplace_back(fThreshAlg->GetPulse(pulseCounter));
 
     // We have to refine this algorithm later
     for (pmtana::pulse_param const& pulse : pulses)
     {
       if (pulse.t_end <= pulse.t_start)
         // Can I call it a logic error?
         throw art::Exception(art::errors::LogicError)
           << "Pulse ends before or at the same time it starts!\n";
 
       std::vector< short >::const_iterator window_start =
               waveform.begin() + static_cast< size_t >(pulse.t_start);
       std::vector< short >::const_iterator window_end   =
               waveform.begin() + static_cast< size_t >(pulse.t_end  );
       mapTickWaveform.emplace(static_cast< size_t >(pulse.t_start),
                               std::vector< short >(window_start, window_end));
       // Don't forget to check that the time output by the (new) algortihm
       // is < fTimeEnd and > fTimeBegin!
 
     }
 
     return mapTickWaveform;
 
   }

double opdet::OpDetDigitizerDUNEDP::SumOfElements ( std::vector< double > aa )

private

Definition at line 962 of file OpDetDigitizerDUNEDP_module.cc.

   {
     double sum=0;
     for (auto& n : aa)    sum+=n;
     return sum;
   }

double opdet::OpDetDigitizerDUNEDP::TickToTime ( size_t tick ) const

private

Definition at line 907 of file OpDetDigitizerDUNEDP_module.cc.

   {
 
 //std::cout << "tick "<< tick<< " " << (tick > fPreTrigger) << " " << (static_cast< double >(tick - fPreTrigger)/fSampleFreq
   //                                                               + fTimeBegin) << " " << (static_cast< double >(fPreTrigger - tick)/fSampleFreq*(-1.0)
 //                                                                 + fTimeBegin)<< std::endl;
     if (tick > fPreTrigger)
       return (static_cast< double >(tick - fPreTrigger)/fSampleFreq
                                                                  + fTimeBegin);
     else
       return (static_cast< double >(fPreTrigger - tick)/fSampleFreq*(-1.0)
                                                                  + fTimeBegin);
 
   }

size_t opdet::OpDetDigitizerDUNEDP::TimeToTick ( double time ) const

private

Definition at line 923 of file OpDetDigitizerDUNEDP_module.cc.

   {
 
     return static_cast< size_t >(std::round((time - fTimeBegin)*fSampleFreq
                                                                + fPreTrigger));
 
   }

std::vector< short > opdet::OpDetDigitizerDUNEDP::VectorOfDoublesToVectorOfShorts ( std::vector< double > const & vectorOfDoubles ) const

private

Definition at line 830 of file OpDetDigitizerDUNEDP_module.cc.

   {
     // Don't bother to round properly, it's faster this way
     return std::vector<short>(vectorOfDoubles.begin(), vectorOfDoubles.end());
 
     /*
     std::vector< short > vectorOfShorts;
     vectorOfShorts.reserve(vectorOfDoubles.size());
     int i=0;
     for (short const& value : vectorOfDoubles)
     {
       vectorOfShorts.emplace_back(static_cast< short >(std::round(value)));
     }
     return vectorOfShorts;
  //   */
   }