33 #include "NovaDAQConventions/DAQConventions.h" 41 #include "Utilities/func/MathUtil.h" 77 double PE2ADC(
double smearedPE,
double gain);
81 std::unique_ptr< std::vector<rawdata::RawDigit> > &rawdigitlist,
168 produces< std::vector<rawdata::RawDigit> >();
169 produces< std::vector<rawdata::RawTrigger> >();
171 std::vector<int> tracePlanes = pset.
get<std::vector<int> >(
"TracePlanes");
172 std::vector<int> traceCells = pset.
get<std::vector<int> >(
"TraceCells");
173 if(traceCells.size() != tracePlanes.size())
175 "FCL configurable TracePlanes and TraceCells must have " <<
176 "same number of elements \n";
178 for(
size_t i = 0;
i < traceCells.size(); ++
i){
179 int plane = tracePlanes[
i];
180 int cell = traceCells[
i];
204 switch(geom->
DetId()){
218 assert(0 &&
"Unknown detector");
245 int runNumber = run.
run();
246 if(runNumber >= 1000000 || runNumber < 10000) {
247 int forceGain = pset.
get<
int>(
"ForceGain");
248 if (forceGain == 0) {
250 <<
"Ideal MC run number detected, but ForceGain == 0. Set ForceGain to the value of gain you want to simulate.\n" 251 << __FILE__ <<
":" << __LINE__ <<
"\n";
266 if(noiseFile.empty() ||
stat(noiseFile.c_str(), &sb) != 0){
269 <<
"noise file " << noiseFile <<
" not found!\n" 270 << __FILE__ <<
":" << __LINE__ <<
"\n";
273 if(gainFile.empty() ||
stat(gainFile.c_str(), &sb) != 0){
276 <<
"gain file " << gainFile <<
" not found!\n" 277 << __FILE__ <<
":" << __LINE__ <<
"\n";
282 mf::LogError(
"ReadoutSimBeginJob") <<
"failed to load cell noise pdf";
287 if (gSystem->AccessPathName(gainFile.c_str())){
289 <<
"Unable to open " << gainFile <<
"\n" 290 << __FILE__ <<
":" << __LINE__ <<
"\n";
293 TFile
f(gainFile.c_str());
294 for (
int ipix = 0; ipix < 32; ++ipix)
296 TString gainName =
"hGains_PIX";
298 TH1D *htemp = (TH1D*)
f.Get(gainName);
300 throw cet::exception(
"FatalRootError") <<
"Unable to get gain histogram from " << gainFile <<
"\n" 301 << __FILE__ <<
":" << __LINE__ <<
"\n";
303 gainName =
"GainPDF_PIX";
305 fGainPDF.push_back( (TH1D*)htemp->Clone(gainName) );
312 std::vector< float > plane_gains;
323 double baselineMean = pset.
get<
double>(
"BaselineMean");
324 double baselineSigma = pset.
get<
double>(
"BaselineSigma");
329 std::vector< double > plane_baselines;
331 double baseline = gauss.
fire(baselineMean, baselineSigma);
332 plane_baselines.push_back(baseline);
348 fADC = tfs->
make<TH1D>(
"ADC",
";ADC", 4095, 0, 4095);
349 fMaxFebADC = tfs->
make<TH1D>(
"MaxFebADC",
";MaxFebADC", 2000, 0, 20000);
350 fNoiseADC = tfs->
make<TH1D>(
"NoiseADC",
";ADC3-ADC0", 150, 0, 300);
363 std::unique_ptr< std::vector<rawdata::RawDigit> > rawdigitlist(
new std::vector<rawdata::RawDigit>);
364 std::unique_ptr< std::vector<rawdata::RawTrigger> > rawtriglist(
new std::vector<rawdata::RawTrigger>);
387 double gainScale =
fGainPDF[pix]->GetMean(1);
397 std::vector<sim::PhotonSignal> hitlist;
398 for(
unsigned int i = 0;
i < hitlistHandle->size(); ++
i){
399 if((*hitlistHandle)[
i].NPhoton() > 0)
400 hitlist.push_back((*hitlistHandle)[
i]);
413 rawtriglist->push_back(rawtrig);
414 evt.
put(std::move(rawtriglist));
419 mf::LogWarning(
"FailedDigitCreation") <<
"failed to create RawDigits, \n" 420 <<
"putting empty RawDigit list into event";
421 evt.
put(std::move(rawdigitlist));
426 for(
unsigned int i = 0;
i < rawdigitlist->size(); ++
i){
427 TString debugOut = cmap->
GetPlane(&rawdigitlist->at(
i));
429 debugOut += cmap->
GetCell(&rawdigitlist->at(
i));
431 for (
int iADC = 0; iADC < (
int)rawdigitlist->at(
i).NADC(); ++ iADC) {
432 debugOut += rawdigitlist->at(
i).ADC(iADC);
438 fADC->Fill(rawdigitlist->at(
i).ADC(3) - rawdigitlist->at(
i).ADC(0));
444 for(
unsigned int i = 0;
i < rawdigitlist->size(); ++
i)
445 rawdigitlist->at(
i).SetMC();
455 if (digitcol->empty())
458 for (
const auto & srcDigit : *digitcol)
461 if (!srcDigit.IsMC())
462 rawdigitlist->push_back(srcDigit);
466 evt.
put(std::move(rawdigitlist));
488 (
const std::vector<sim::PhotonSignal>& hitlist,
489 std::unique_ptr< std::vector<rawdata::RawDigit> > &rawdigitlist,
543 std::map< std::pair<int, int>, std::list< std::pair<double, double> > > ADCPulseMap;
545 std::map< std::pair<int, int>, std::list< std::pair<double, double> > > SaggedPulseMap;
547 std::vector<bool> channelHasHits[geom->
NPlanes()];
550 channelHasHits[
plane].push_back(
false);
555 const int plane = phot.Plane();
556 const int cell = phot.Cell();
558 std::pair<int, int> hitCoord(plane, cell);
559 std::pair<int, int> febCoord(plane, cell/32);
561 int nPE = phot.NPhoton();
564 double time = phot.TimeMean();
573 else ADCPulseMap[hitCoord].push_back(
std::make_pair(time, ADC) );
576 for (
auto&&
it: ADCPulseMap)
it.second.sort();
577 for (
auto&&
it: SaggedPulseMap)
it.second.sort();
580 std::vector<bool> febHasFlashHit[geom->
NPlanes()];
584 double meanTime(0.0);
587 for (
auto it = SaggedPulses.begin();
it != SaggedPulses.end(); ++
it) {
588 double currTime =
it->first;
589 double currADC =
fabs(
it->second);
594 else if (
std::fabs(meanTime - currTime) < 15.0 ) {
595 meanTime = meanTime*ADC + currTime*currADC;
600 if (ADC > maxADC) maxADC = ADC;
605 if (ADC > maxADC) maxADC = ADC;
609 else febHasFlashHit[
plane].push_back(
false);
628 std::vector<std::vector<int>> phases;
629 phases.resize(geom->
NPlanes());
630 for(
unsigned int i = 0;
i < geom->
NPlanes(); ++
i){
647 const int phase = phases[
plane][
cell/32];
658 double noiseFactor = 1.0;
690 for (
int iClocktick=firstDigitizationClocktick;
694 ASIC_Output[iClocktick] =
int(ASIC_Output[iClocktick]);
695 if ( ASIC_Output[iClocktick] < 0 ) ASIC_Output[iClocktick] = 0 ;
696 else if ( ASIC_Output[iClocktick] > (1<<12)-1 ) ASIC_Output[iClocktick] = (1<<12)-1;
702 if(
plane != (
unsigned int)tracePlaneCell.first)
continue;
704 if(
cell/32 != (
unsigned int)tracePlaneCell.second/32)
continue;
724 std::vector<rawdata::RawDigit> local_digitlist =
726 LOG_DEBUG(
"ReadoutSim") <<
" == size of digilist == " 727 << local_digitlist.size();
731 for (
unsigned int iD=0; iD<local_digitlist.size(); iD++) {
734 if ( ! digit.
NADC() ) {
736 <<
"ReadoutSim::CreateRawDigits() -- FPGA algorithm returned " 737 <<
"a digit with no ADC values. This should never happen. Fatal!\n" 738 << __FILE__ <<
":" << __LINE__ <<
"\n";
747 <<
" " << digit.
TDC() <<
" " << digit.
ADC(0);
749 rawdigitlist->push_back(digit);
757 double mean_arrival = 0;
759 mean_arrival = 1.0/
mean;
762 for (Long64_t iTick = firstDigitizationClocktick;
768 if (iTick < fParams->fNumClockticksInSpill) {
770 double adc0, adc1, adc2;
784 else digit.
SetADC(0,uint16_t(adc0));
787 digit.
SetADC(0,uint16_t(baseline));
788 digit.
SetADC(1,uint16_t(adc0 + baseline));
789 digit.
SetADC(2,uint16_t(adc1 + baseline));
790 digit.
SetADC(3,uint16_t(adc2 + baseline));
797 if (digit.
ADC() >= thresh) {
799 rawdigitlist->push_back(digit);
840 <<
"ReadoutSim::GetASICScaleSpread() -- nPE==0. This should never happen.\n" 841 << __FILE__ <<
":" << __LINE__ <<
"\n";
849 const double mu = -
log(1+var)/2;
853 return nPE*
exp(mu+sigma*gauss.
fire(0, 1));
866 if (gSystem->AccessPathName(fname)) {
867 cet::exception(
"FileNotOpen") <<
"ReadoutSim: Unable to open " << fname
868 <<
" Cannot get noise PDFs!";
874 TH3D *htemp = (TH3D*)f.Get(
"noise/h3_adc0adc1adc2");
876 std::cerr <<
"ReadoutSim: Histogram h3_adc0adc1adc2 not found in " << fname <<
std::endl;
924 static const unsigned long long ticksPerYear = (
unsigned long long)(64000000) * 31556926;
925 static const unsigned long long ticksPerDay = (
unsigned long long)(64000000) * 86400;
926 static const unsigned long long ticksPerEvent = (
unsigned long long)(64000000) / 10;
929 evt.
run() * ticksPerYear +
930 evt.
subRun() * ticksPerDay +
931 evt.
event() * ticksPerEvent;
944 std::cout <<
"Saving a trace for plane: " << plane
947 std::vector<double> traceTimes;
948 std::vector<double>
trace;
950 for(
int iClocktick = offset;
956 trace.push_back(ASIC_Output[iClocktick]);
961 TGraph* traceGraph = tfs->
make<TGraph>(traceTimes.size(),
962 &traceTimes[0], &trace[0]);
963 std::stringstream traceName;
964 traceName<<name<<
"_plane" <<plane<<
"_feb"<<cell/32<<
"_cell"<<
cell;
965 traceGraph->SetName(traceName.str().c_str());
966 traceGraph->SetTitle(
";Time [ns]; ADC");
double fAPDExcessNoiseFactor
APD's "excess noise factor" (see .cxx for more)
void SortByPlaneAndCell(std::vector< sim::PhotonSignal > &c)
#define LOG_DEBUG(stream)
int fNumClockticksPerDigitization
How many ticks between each ADC.
SubRunNumber_t subRun() const
back track the reconstruction to the simulation
double DrawSmearedPE(int nPE, CLHEP::RandFlat *flat)
virtual std::vector< double > CreateTrace(std::list< std::pair< double, double > > &ADCPulses, std::list< std::pair< double, double > > &SaggedPulses)=0
virtual std::vector< rawdata::RawDigit > ASICToDigits(double *ASIC_Output, int offset)=0
int32_t TDC() const
The time of the last baseline sample.
fvar< T > fabs(const fvar< T > &x)
std::vector< TH1 * > fGainPDF
void beginRun(art::Run &run)
double GetNominalThreshold()
void SetTDC(int32_t iTDC)
ReadoutSim(const fhicl::ParameterSet &pset)
rsim::IPulseShaper * fPulseShaper
PulseShaper object for generating traces.
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
TH3D * fEmptyCellNoisePDF
for fast generation of empty-cell noise hits
std::string fEmptyCellNoiseFile_NDOS
File holding NDOS noise PDF.
bool fGeneratePhysicsCellNoise
void SetBaseline(double Baseline)
unsigned int Ncells() const
Number of cells in this plane.
base_engine_t & createEngine(seed_t seed)
void SetChannel(uint32_t iChan)
const daqchannelmap::DAQChannelMap * Map() const
::xsd::cxx::tree::exception< char > exception
IFPGAAlgorithm * CreateFPGAAlgorithm(const fhicl::ParameterSet &pset)
Create and configure the correct algorithm based on pset.
void PopulateRawTrig(rawdata::RawTrigger &rawtrig, art::Event const &evt)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
std::string fGainFile
File containing the distribution of gains on the detector.
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
std::string find_file(std::string const &filename) const
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
std::string fEmptyCellNoiseFile_FD
File holding FD noise PDF.
rsim::ExcessNoiseMaker * fExcessNoiseMaker
ExcessNoiseMaker object for modeling.
std::string fPhotonModuleLabel
photons made in the module with this label
std::string fEmptyCellNoiseFile_TB
File holding TB noise PDF.
void SetVersion(uint8_t v)
rsim::IFPGAAlgorithm * fFPGAAlgo
ProductID put(std::unique_ptr< PROD > &&product)
double fScaleEmptyCellNoiseRate
scale factor for empty cell noise
Common configuration params for SimpleReadout, FPGAAlgorithms, NoiseMaker.
int LoadEmptyCellNoisePDF(const char *fname)
void SetADC(uint32_t i, int16_t iADC)
base_engine_t & getEngine() const
double GetNoise(CLHEP::RandGaussQ *gauss)
double fFEBFlashThreshold
void SetPlaneAndCell(int plane, int cell)
Far Detector at Ash River, MN.
fhicl::ParameterSet fPSetND
std::string fCopyNoTruthHitsLabel
hits with no truth (i.e., overlay) from this label will be copied into the digit vector after simulat...
bool fUseNewEmptyCellNoise
For backwards compatibility with the old noise model until the 2nd analysis is over.
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
fvar< T > exp(const fvar< T > &x)
Interface implemented by all FPGA algorithms.
Near Detector in the NuMI cavern.
std::vector< std::vector< float > > fGainMap
unsigned short GetPlane(const rawdata::RawDigit *dig)
double fGain
APD gain (electrons per photoelectron)
void SetDaqChannel(uint32_t iChan)
bool CreateRawDigits(const std::vector< sim::PhotonSignal > &hitlist, std::unique_ptr< std::vector< rawdata::RawDigit > > &rawdigitlist, rawdata::RawTrigger rawtrig)
void produce(art::Event &evt)
fhicl::ParameterSet fPSetTB
double ThresholdADC() const
EventNumber_t event() const
std::vector< std::vector< double > > fBaselineMap
fhicl::ParameterSet fPSetNDOS
static constexpr Double_t gauss
rsim::NoiseMaker * fNoiseMaker
NoiseMaker object for making correlated noise.
void Rebuild(const art::Event &evt)
double fADCMaxPE
maximum signal size that the ADC can handle without saturating (in PE)
std::vector< std::pair< int, int > > fTracePlaneCells
unsigned int GetRandomNumberSeed()
pixel_t getPixel(dchan daqchan) const
Decode the pixel id from a dchan.
double fEmptyCellNoisePDFIntegral
bool fVaryNoiseByThreshold
fhicl::ParameterSet fPSetFD
T * make(ARGS...args) const
std::string fEmptyCellNoiseFile_ND
File holding ND noise PDF.
int16_t ADC(uint32_t i) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
unsigned long long fTriggerTimingMarker_TimeStart
fvar< T > floor(const fvar< T > &x)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Can be used as either a member holding configurations, or a mix-in.
T trace(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
int fNumClockticksInSpill
number of clockticks in a spill
double GetSmearedPE(int nPE)
assert(nhit_max >=nhit_nbins)
unsigned short GetCell(const rawdata::RawDigit *dig)
void GetRandom(TH1 *histo, double &x, double &y, double &z)
uint32_t fTriggerRange_TriggerLength
dchan encodeDChan(int detID, diblock_t diblock, dcm_id_t dcm, feb_t feb, pixel_t pixel) const
CommonParameters * fParams
unsigned int NPlanes() const
double PE2ADC(double smearedPE, double gain)
uint32_t dchan
< DAQ Channel Map Package
std::string fEmptyCellNoiseFile
File containing the empty cell noise templates.
Encapsulate the geometry of one entire detector (near, far, ndos)
double fASICBaselineInADCCounts
ASIC baseline output level, in ADC counts.
double fClocktick
digitization clock period, ns
void DrawASICOutput(double *ASIC_Output, int offset, int stride, int plane, int cell, std::string name="trace_")