NoiseCorrelation_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: NoiseCorrelation
3 // Module Type: analyser
4 // File: NoiseCorrelation_module.cc
5 // Author: Mike Wallbank (m.wallbank@sheffield.ac.uk), Feb 2016
6 //
7 // Looks at the correlation between waveforms on all combinations
8 // of channels.
9 // Runs over an artdaq-formatted file and produces a tree or histogram,
10 // filled for each channel combination and the respective correlation.
11 ////////////////////////////////////////////////////////////////////////
12 
13 // framework
20 #include "art_root_io/TFileService.h"
21 #include "art_root_io/TFileDirectory.h"
23 #include "fhiclcpp/ParameterSet.h"
25 
26 // lbne-artdaq
27 #include "lbne-raw-data/Overlays/TpcMilliSliceFragment.hh"
28 #include "artdaq-core/Data/Fragment.hh"
29 
30 // larsoft
32 #include "lardataobj/RawData/raw.h"
34 #include "tpcFragmentToRawDigits.h"
41 #include "cetlib/getenv.h"
42 
43 // c++
44 #include <memory>
45 #include <array>
46 #include <iostream>
47 #include <fstream>
48 #include <sstream>
49 
50 // ROOT
51 #include "TMath.h"
52 #include "TTree.h"
53 #include "TH2F.h"
54 #include "TStyle.h"
55 
56 namespace DAQToOffline {
57  class NoiseCorrelation;
58 }
59 
61 public:
62 
63  explicit NoiseCorrelation(fhicl::ParameterSet const & pset);
64  // The destructor generated by the compiler is fine for classes
65  // without bare pointers or other resource use.
66 
67  // Plugins should not be copied or assigned.
68  NoiseCorrelation(NoiseCorrelation const &) = delete;
70  NoiseCorrelation & operator = (NoiseCorrelation const &) = delete;
72 
73  void analyze(art::Event const& evt) override;
74  void reconfigure(const fhicl::ParameterSet &pset);
75 
76 private:
77 
80 
81  // Correlation tree
84  float fCorrelation;
85  int fChannel1;
86  int fChannel2;
87 
88  const lariov::DetPedestalProvider& fPedestalRetrievalAlg = *(lar::providerFrom<lariov::DetPedestalService>());
89 
90 };
91 
93 
94  this->reconfigure(pset);
95  gStyle->SetOptStat(0);
96 
98  fCorrelationTree = tfs->make<TTree>("Correlations","Correlations");
99  fCorrelationHist = tfs->make<TH2F>("Correlation","Correlation;Online channel number;Online channel number;",2048,0,2048,2048,0,2048);
100  fCorrelationTree->Branch("Channel1",&fChannel1);
101  fCorrelationTree->Branch("Channel2",&fChannel2);
102  fCorrelationTree->Branch("Correlation",&fCorrelation);
103 
104 }
105 
107  fFragType = pset.get<std::string>("FragType");
108  fRawDataLabel = pset.get<std::string>("RawDataLabel");
109 }
110 
112 
113  art::Handle<artdaq::Fragments> rawFragments;
114  evt.getByLabel(fRawDataLabel, fFragType, rawFragments);
115 
116  art::EventNumber_t eventNumber = evt.event();
117 
118  // Check that the data are valid
119  if(!rawFragments.isValid()){
120  std::cerr << "Run: " << evt.run()
121  << ", SubRun: " << evt.subRun()
122  << ", Event: " << eventNumber
123  << " is NOT VALID" << std::endl;
124  throw cet::exception("rawFragments NOT VALID");
125  }
126 
127  // Create a map containing (fragmentID, fragIndex) for the event, will be used to check if each channel is present
128  std::map < unsigned int, unsigned int > mapFragID;
129  for(size_t fragIndex = 0; fragIndex < rawFragments->size(); fragIndex++){
130  const artdaq::Fragment &singleFragment = (*rawFragments)[fragIndex];
131  unsigned int fragmentID = singleFragment.fragmentID();
132  mapFragID.insert(std::pair<unsigned int, unsigned int>(fragmentID,fragIndex));
133  }
134 
135  // Create a 2D vector to save the correlations
136  std::vector<std::vector<float> > correlationArray(2048,std::vector<float>(2048,0));
137 
139  size_t numChannels = geometry->Nchannels();
140 
141  // Get the ADC vector for each channel
142  std::vector<std::vector<short> > adcVectors;
143 
144  for (size_t channel = 0; channel < numChannels; ++channel) {
145 
146  // Follow the steps written by J Davies in TpcDAQToOffline
147  unsigned int fragmentID = UnpackFragment::getFragIDForChan(channel);
148  unsigned int sample = UnpackFragment::getNanoSliceSampleForChan(channel);
149  std::vector<short> adcvec;
150  if (mapFragID.find(fragmentID) != mapFragID.end()) {
151  unsigned int fragIndex = mapFragID[fragmentID];
152  const artdaq::Fragment &singleFragment = (*rawFragments)[fragIndex];
153  lbne::TpcMilliSliceFragment millisliceFragment(singleFragment);
154  auto numMicroSlices = millisliceFragment.microSliceCount();
155  for(unsigned int i_micro = 0; i_micro < numMicroSlices; i_micro++) {
156  std::unique_ptr <const lbne::TpcMicroSlice> microSlice = millisliceFragment.microSlice(i_micro);
157  auto numNanoSlices = microSlice->nanoSliceCount();
158  for(uint32_t i_nano = 0; i_nano < numNanoSlices; i_nano++){
160  bool success = microSlice->nanosliceSampleValue(i_nano, sample, val);
161  if (success) {
162  float const pedestal = fPedestalRetrievalAlg.PedMean(channel);
163  short adc = short(val - pedestal);
164  if ((adc & 0x3F) == 0x0 || (adc & 0x3F) == 0x3F)
165  adcvec.push_back(short(-999));
166  else
167  adcvec.push_back(adc);
168  }
169  }
170  }
171  }
172 
173  adcVectors.push_back(adcvec);
174 
175  }
176 
177  // Look at all possible channel combinations
178  for (size_t chan1 = 0; chan1 < numChannels; ++chan1) {
179  for (size_t chan2 = chan1; chan2 < numChannels; ++chan2) {
180 
181  fChannel1 = chan1;
182  fChannel2 = chan2;
183  fCorrelation = -999;
184 
185  std::vector<short> adcvec1 = adcVectors[chan1];
186  std::vector<short> adcvec2 = adcVectors[chan2];
187 
188  if (fChannel1 % 50 == 0 and fChannel2 % 1000 == 0)
189  std::cout << "Looking at correlations between channel " << fChannel1 << " and " << fChannel2 << std::endl;
190 
191  // Now we have the ADC vectors, we can look at the correlations
192  // Use the Pearson Correlation Coefficient
193 
194  size_t n1 = adcvec1.size(), n2 = adcvec2.size();
195  if (n1 != n2)
196  continue;
197 
198  float sumxy = 0, sumx = 0, sumy = 0, sumx2 = 0, sumy2 = 0;
199  for (size_t tick = 0; tick < n1; ++tick) {
200  short adc1 = adcvec1[tick], adc2 = adcvec2[tick];
201  sumx += adc1;
202  sumy += adc2;
203  sumxy += adc1*adc2;
204  sumx2 += adc1*adc1;
205  sumy2 += adc2*adc2;
206  }
207 
208  float denom = ( (n1*sumx2) - (sumx*sumx) ) * ( (n1*sumy2) - (sumy*sumy) );
209  if (n1 > 0 and denom > 0) {
210  fCorrelation = ( (n1*sumxy) - (sumx*sumy) ) / ( TMath::Sqrt( denom ) );
211  correlationArray[fChannel1][fChannel2] = fCorrelation;
212  correlationArray[fChannel2][fChannel1] = fCorrelation;
213  }
214  //fCorrelationTree->Fill();
215 
216  }
217  } // channel loops
218 
219  for (unsigned int channel1 = 0; channel1 < correlationArray.size(); ++channel1)
220  for (unsigned int channel2 = 0; channel2 < correlationArray.size(); ++channel2)
221  fCorrelationHist->SetBinContent(channel1, channel2, correlationArray[channel1][channel2]);
222  fCorrelationHist->GetZaxis()->SetRangeUser(-1,1);
223 
224  return;
225 
226 }
227 
EventNumber_t event() const
Definition: DataViewImpl.cc:96
std::string string
Definition: nybbler.cc:12
unsigned int getNanoSliceSampleForChan(unsigned int channel)
Definition of basic raw digits.
unsigned short uint16_t
Definition: stdint.h:125
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:27
NoiseCorrelation & operator=(NoiseCorrelation const &)=delete
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
unsigned int getFragIDForChan(unsigned int channel)
bool isValid() const
Definition: Handle.h:183
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:69
unsigned int uint32_t
Definition: stdint.h:126
void analyze(art::Event const &evt) override
Collect all the RawData header files together.
T get(std::string const &key) const
Definition: ParameterSet.h:231
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:74
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:89
RunNumber_t run() const
Definition: DataViewImpl.cc:82
static int max(int a, int b)
NoiseCorrelation(fhicl::ParameterSet const &pset)
void reconfigure(const fhicl::ParameterSet &pset)
const lariov::DetPedestalProvider & fPedestalRetrievalAlg
IDNumber_t< Level::Event > EventNumber_t
Definition: IDNumber.h:118
Interface for experiment-specific channel quality info provider.
virtual float PedMean(raw::ChannelID_t ch) const =0
Retrieve pedestal information.
TCEvent evt
Definition: DataStructs.cxx:7
Interface for experiment-specific service for channel quality info.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)