MakeNoiseSpectrumFile_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief TODO
3 /// \author Christopher Backhouse
4 ////////////////////////////////////////////////////////////////////////
5 
8 
9 // Framework includes
16 
17 #include "fhiclcpp/ParameterSet.h"
18 
20 
21 // NOvA includes
22 #include "Simulation/Simulation.h"
23 
24 #include "TH3.h"
25 #include "TMatrixD.h"
26 #include "TVector.h"
27 
28 namespace rsim
29 {
30  const int kCovMatSize = 6;
31 
33  {
34  public:
35  explicit MakeNoiseSpectrumFile(fhicl::ParameterSet const &pset);
36  virtual ~MakeNoiseSpectrumFile();
37 
38  void analyze(const art::Event& evt);
39  void beginJob();
40  void endJob();
41 
42  protected:
43  int fN; ///< How many spills have we simulated?
44 
45  // Information for calculating the covariance matrix of the noise.
46  // Need to use TMatrixD because default TMatrix overflows.
51 
54 
55  TH3* fNoiseMap;
56  };
57 
58  //----------------------------------------------------------------------
60  EDAnalyzer(pset),
61  CommonParameters(pset),
62  fN(0),
63  fCovMat(kCovMatSize, kCovMatSize), fCovMatStats(kCovMatSize, kCovMatSize),
64  fMeanVec(kCovMatSize), fMeanVecStats(kCovMatSize),
66  fNoiseMaker(pset)
67  {
68  // get the random number seed, use a random default if not specified
69  // in the configuration file.
70  unsigned int seed = pset.get<unsigned int>("Seed",
72 
73  createEngine(seed);
74  }
75 
76  //----------------------------------------------------------------------
78  {
79  delete fFPGAAlgo;
80  }
81 
82  //----------------------------------------------------------------------
84  {
86 
87  fNoiseMap = tfs->make<TH3D>("h3_adc0adc1adc2", ";ADC0;ADC1;ADC2",
88  80, 0, 80,
89  80, 0, 80,
90  80, 0, 80);
91  }
92 
93  //----------------------------------------------------------------------
95  {
96  ++fN;
97 
99  std::cout.flush();
101 
102  const int maxADC = (1<<12)-1;
103  double ASIC_Output[fNumClockticksInSpill];
104 
105  for(int n = 0; n < fNumClockticksInSpill;
107  ASIC_Output[n] = int(fASICBaselineInADCCounts +
109 
110  if(ASIC_Output[n] < 0) ASIC_Output[n] = 0;
111  if(ASIC_Output[n] > maxADC) ASIC_Output[n] = maxADC;
112  }
113 
114 
115  // Fill covariance matrix information
117  for(int i = 0; i < kCovMatSize; ++i){
118  const int iidx = n+fNumClockticksPerDigitization*i;
119  if(iidx >= fNumClockticksInSpill) continue;
120 
121  fMeanVec(i) += ASIC_Output[iidx];
122  ++fMeanVecStats(i);
123 
124  for(int j = 0; j < kCovMatSize; ++j){
125  const int jidx = n+fNumClockticksPerDigitization*j;
126  if(jidx >= fNumClockticksInSpill) continue;
127 
128  fCovMat(i, j) += ASIC_Output[iidx]*ASIC_Output[jidx];
129  ++fCovMatStats(i, j);
130  } // end for j
131  } // end for i
132  } // end for n
133 
134 
135  const std::vector<rawdata::RawDigit> digits =
136  fFPGAAlgo->ASICToDigits(ASIC_Output, 0);
137 
138  for(unsigned int i = 0; i < digits.size(); ++i){
139  const rawdata::RawDigit& digit = digits[i];
140 
141  fNoiseMap->Fill(digit.ADC(0), digit.ADC(1), digit.ADC(2));
142  }
143  }
144 
145  //----------------------------------------------------------------------
147  {
148  // Normalization of noise map is probability /cell /ns
149  // Warning, this can overflow if performed with ints.
150  const double ns = long(fN)*long(fNumClockticksInSpill)*fClocktick;
151  fNoiseMap->Scale(1/ns);
152 
153  // Divide by number of samples to get means
154  for(int i = 0; i < kCovMatSize; ++i){
155  fMeanVec(i) /= fMeanVecStats(i);
156  for(int j = 0; j < kCovMatSize; ++j){
157  fCovMat(i, j) /= fCovMatStats(i, j);
158  }
159  }
160 
161  // Definition of covariance: <ij>-<i><j>
162  for(int i = 0; i < kCovMatSize; ++i){
163  for(int j = 0; j < kCovMatSize; ++j){
164  fCovMat(i, j) -= fMeanVec(i)*fMeanVec(j);
165  }
166  }
167 
168  std::cout << "Covariance matrix: " << std::endl;
169  fCovMat.Print();
170  std::cout << "And its inverse: " << std::endl;
171  fCovMat.Invert();
172  fCovMat.Print();
173 
174  /*
175  // Output format useful for copy-pasting into code
176  for(int i = 0; i < kCovMatSize; ++i){
177  for(int j = 0; j < kCovMatSize; ++j){
178  std::cout << "mat(" << i << ", " << j << ") = "
179  << fCovMat(i, j) << "; ";
180  }
181  std::cout << std::endl;
182  }
183  */
184  }
185 
187 
188 } // namespace
int fNumClockticksPerDigitization
How many ticks between each ADC.
virtual std::vector< rawdata::RawDigit > ASICToDigits(double *ASIC_Output, int offset)=0
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
base_engine_t & createEngine(seed_t seed)
IFPGAAlgorithm * CreateFPGAAlgorithm(const fhicl::ParameterSet &pset)
Create and configure the correct algorithm based on pset.
DEFINE_ART_MODULE(TestTMapFile)
int fN
How many spills have we simulated?
Common configuration params for SimpleReadout, FPGAAlgorithms, NoiseMaker.
base_engine_t & getEngine() const
double GetNoise(CLHEP::RandGaussQ *gauss)
Definition: NoiseMaker.h:32
unsigned int seed
Definition: runWimpSim.h:102
T get(std::string const &key) const
Definition: ParameterSet.h:231
Interface implemented by all FPGA algorithms.
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
static constexpr Double_t gauss
Definition: Munits.h:360
static const double ns
Module that plots metrics from reconstructed cosmic ray data.
OStream cout
Definition: OStream.cxx:6
void analyze(const art::Event &evt)
double fADCMaxPE
maximum signal size that the ADC can handle without saturating (in PE)
unsigned int GetRandomNumberSeed()
Definition: Simulation.cxx:13
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
MakeNoiseSpectrumFile(fhicl::ParameterSet const &pset)
Can be used as either a member holding configurations, or a mix-in.
int fNumClockticksInSpill
number of clockticks in a spill
double fASICBaselineInADCCounts
ASIC baseline output level, in ADC counts.
double fClocktick
digitization clock period, ns