LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
CalWireAna_module.cc
Go to the documentation of this file.
1 //
3 // CalWireAna class designed to make histos.
4 //
5 //
6 //
8 
10 
11 // C++ includes
12 #include <algorithm>
13 #include <string>
14 
15 // Framework includes
20 #include "art_root_io/TFileService.h"
23 #include "fhiclcpp/ParameterSet.h"
25 
26 // LArSoft includes
31 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
32 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
33 
34 // ROOT includes
35 #include "TComplex.h"
36 #include "TH1.h"
37 #include "TH2.h"
38 #include "TMath.h"
39 
40 namespace caldata {
41 
43  class CalWireAna : public art::EDAnalyzer {
44 
45  public:
46  explicit CalWireAna(fhicl::ParameterSet const& pset);
47 
48  private:
50  void analyze(const art::Event& evt) override;
51  void beginJob() override;
52 
53  std::string fCalWireModuleLabel;
54  std::string fDetSimModuleLabel; //< name of module that produced the digits
55 
56  TH1F* fDiffsR;
57  TH1F* fDiffsW;
58  TH1F* fDiffsRW;
59  TH1F* fDiffsRWgph;
60  TH1F* fMin;
61  TH1F* fMax;
62  TH1F* fIR;
63  TH1F* fCR;
64  TH1F* fIW;
65  TH1F* fCW;
66  TH1F* fRawIndPeak; //Raw Induction Peak Values
67  TH1F* fRawColPeak; //Raw Collection Peak Values
68  TH1F* fCalIndPeak; //Calibrated Induction Peak Values
69  TH1F* fCalColPeak; //Calibrated Collection Peak Values
70  TH1F* fNoiseHist; //Noise Frequency Spectrum
71  TH1F* fNoiseRMS; //Noise RMS values
72  TH2F* fWireSig;
73  TH2F* fRawSig;
74  TH2F*
76  TH2F*
78  TH2F* fDiffsRWvsR;
80  TH2F* fWindow;
81 
82  }; // class CalWireAna
83 
84 } // End caldata namespace.
85 
86 namespace caldata {
87 
88  //-------------------------------------------------
90  : EDAnalyzer(pset)
91  , fCalWireModuleLabel(pset.get<std::string>("CalWireModuleLabel"))
92  , fDetSimModuleLabel(pset.get<std::string>("DetSimModuleLabel"))
93  {}
94 
95  //-------------------------------------------------
97  {
98  // get access to the TFile service
100 
101  fDiffsR = tfs->make<TH1F>("One timestamp diffs in RawDigits", ";#Delta ADC;", 40, -19.5, 20.5);
102  fDiffsW = tfs->make<TH1F>("One timestamp diffs in Wires", ";#Delta ADC;", 20, -9.5, 10.5);
103  fDiffsRW = tfs->make<TH1F>("Same timestamp diffs in RD-Wires", ";#Delta ADC;", 20, -9.5, 10.5);
104  fDiffsRWvsR = tfs->make<TH2F>(
105  "Same timestamp diffs in RD-Wires vs R", ";#Delta ADC;", 481, -0.5, 480.5, 20, -9.5, 10.5);
106  fDiffsRWgph =
107  tfs->make<TH1F>("Same timestamp diffs in RD-Wires gph", ";#Delta ADC;", 20, -9.5, 10.5);
108  fDiffsRWvsRgph = tfs->make<TH2F>("Same timestamp diffs in RD-Wires vs R gph",
109  ";#Delta ADC;",
110  481,
111  -0.5,
112  480.5,
113  20,
114  -9.5,
115  10.5);
116  fRawSig =
117  tfs->make<TH2F>("One event, one channel Raw", "timestamp", 481, -0.5, 480.5, 21, -0.5, 20.5);
118  fWireSig =
119  tfs->make<TH2F>("One event, one channel Wire", "timestamp", 481, -0.5, 480.5, 21, -0.5, 20.5);
120 
121  fRD_WireMeanDiff2D = tfs->make<TH2F>(
122  "Mean (Raw-CALD)-over-Raw in Window gph", "Wire number", 481, -0.05, 480.5, 40, -1., 1.);
123  fRD_WireRMSDiff2D = tfs->make<TH2F>(
124  "RMS (Raw-CALD)-over-Raw in Window gph", "Wire number", 481, -0.05, 480.5, 10, 0., 2.);
125 
126  fWindow = tfs->make<TH2F>("tmax-tmin vs indMax", "ticks", 200, 0, 2000, 20, -2.5, 60.5);
127  fMin = tfs->make<TH1F>("Value of min", "ticks", 21, -20.5, 0.5);
128  fMax = tfs->make<TH1F>("Value of max", "ticks", 21, 0.5, 20.5);
129  fRawIndPeak = tfs->make<TH1F>("indPeakRaw", ";Induction Peaks Raw;", 40, 5, 45);
130  fRawColPeak = tfs->make<TH1F>("colPeakRaw", ";Collection Peaks Raw;", 40, 5, 45);
131  fCalIndPeak = tfs->make<TH1F>("indPeakCal", ";Induction Peaks Calibrated;", 40, 5, 45);
132  fCalColPeak = tfs->make<TH1F>("colPeakCal", ";Collection Peaks Calibrated;", 40, 5, 45);
133 
134  fIR = tfs->make<TH1F>("Raw Ind signal", "time ticks", 4096, 0.0, 4096.);
135  fCR = tfs->make<TH1F>("Raw Coll signal", "time ticks", 4096, 0.0, 4096.);
136  fIW = tfs->make<TH1F>("Wire Ind signal", "time ticks", 4096, 0.0, 4096.);
137  fCW = tfs->make<TH1F>("Wire Coll signal", "time ticks", 4096, 0.0, 4096.);
138  fNoiseHist = tfs->make<TH1F>("Noise Histogram", "FFT Bins", 2049, 0, 2049);
139  fNoiseRMS = tfs->make<TH1F>("Noise RMS", "RMS", 25, 0, 2.0);
140  }
141 
142  //-------------------------------------------------
144  {
145  // loop over the raw digits and get the adc vector for each, then compress it and uncompress it
146 
147  lariov::ChannelStatusProvider const& channelStatus =
150  evt.getByLabel(fDetSimModuleLabel, rdHandle);
152  evt.getByLabel(fCalWireModuleLabel, wHandle);
153 
154  // get the raw::RawDigit associated by fCalWireModuleLabel to wires in wHandle;
155  // RawDigitsFromWire.at(index) will be a art::Ptr<raw::RawDigit>
156  art::FindOneP<raw::RawDigit> RawDigitsFromWire(wHandle, evt, fCalWireModuleLabel);
157 
159  for (unsigned int i = 0; i < wHandle->size(); ++i) {
160  wvec.emplace_back(wHandle, i);
161  }
163  for (unsigned int i = 0; i < rdHandle->size(); ++i) {
164  rdvec.emplace_back(rdHandle, i);
165  }
166 
168  double pedestal = rdvec[0]->GetPedestal();
169  double threshold = 9.0;
170  double signalSize = rdvec[0]->Samples();
171  uint32_t indChan0 = 64;
172  uint32_t indChan1 = 110;
173  uint32_t colChan0 = 312;
174  uint32_t colChan1 = 354;
175  std::vector<double> ir(fft->FFTSize()), iw(fft->FFTSize()), cr(fft->FFTSize()),
176  cw(fft->FFTSize());
177  ir[signalSize] = iw[signalSize] = cr[signalSize] = cw[signalSize] = 1.0;
179  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
180  for (unsigned int rd = 0; rd < rdvec.size(); ++rd) {
181  // Find corresponding wire.
182  std::vector<double> signal(fft->FFTSize());
183  for (unsigned int wd = 0; wd < wvec.size(); ++wd) {
184  if (RawDigitsFromWire.at(wd) == rdvec[rd]) {
185  std::vector<float> wirSig = wvec[wd]->Signal();
186  if (wirSig.size() > signal.size()) {
187  MF_LOG_DEBUG("CalWireAna")
188  << "Incompatible vector size " << wirSig.size() << " " << signal.size();
189  return;
190  }
191  for (unsigned int ii = 0; ii < wirSig.size(); ++ii) {
192  signal[ii] = wirSig[ii];
193  }
194  break;
195  }
196  if (wd == (wvec.size() - 1)) {
197  MF_LOG_DEBUG("CalWireAna")
198  << "caldata::CalWireAna:Big problem! No matching Wire for RawDigit. Bailing." << rd;
199  return;
200  }
201  }
202 
203  std::vector<double> adc(fft->FFTSize());
204 
205  for (unsigned int t = 1; t < rdvec[rd]->Samples(); ++t) {
206  fDiffsR->Fill(rdvec[rd]->ADC(t) - rdvec[rd]->ADC(t - 1));
207  adc[t - 1] = rdvec[rd]->ADC(t - 1);
208  fRawSig->Fill(rd, rdvec[rd]->ADC(t));
209  }
210  //get the last one for the adc vector
211  adc[rdvec[rd]->Samples() - 1] = rdvec[rd]->ADC(rdvec[rd]->Samples() - 1);
212  if (!channelStatus.IsBad(rdvec[rd]->Channel()) &&
213  (*max_element(adc.begin(), adc.end()) < pedestal + threshold &&
214  *min_element(adc.begin(), adc.end()) > pedestal - threshold)) {
215  double sum = 0;
216  for (int i = 0; i < signalSize; i++)
217  sum += pow(adc[i] - pedestal, 2.0);
218  fNoiseRMS->Fill(TMath::Sqrt(sum / (double)signalSize));
219  std::vector<double> temp(fft->FFTSize());
220  std::vector<TComplex> fTemp(fft->FFTSize() / 2 + 1);
221  for (int i = 0; i < signalSize; i++)
222  temp[i] = (adc[i] - pedestal) * sin(TMath::Pi() * (double)i / signalSize);
223  fft->DoFFT(temp, fTemp);
224  for (int i = 0; i < fft->FFTSize() / 2 + 1; i++)
225  fNoiseHist->Fill(i, fTemp[i].Rho());
226  }
227  if (wireReadoutGeom.SignalType(rdvec[rd]->Channel()) == geo::kInduction &&
228  rdvec[rd]->Channel() > indChan0 && rdvec[rd]->Channel() < indChan1) {
229  fft->AlignedSum(ir, adc);
230  fft->AlignedSum(iw, signal);
231  }
232  if (wireReadoutGeom.SignalType(rdvec[rd]->Channel()) == geo::kCollection &&
233  rdvec[rd]->Channel() > colChan0 && rdvec[rd]->Channel() < colChan1) {
234  fft->AlignedSum(cr, adc);
235  fft->AlignedSum(cw, signal);
236  }
237  if (wireReadoutGeom.SignalType(rdvec[rd]->Channel()) == geo::kInduction) {
238  if (*max_element(adc.begin(), adc.end()) > pedestal + threshold)
239  fRawIndPeak->Fill(*max_element(adc.begin(), adc.end()));
240  if (*max_element(signal.begin(), signal.end()) > pedestal + threshold)
241  fCalIndPeak->Fill(*max_element(signal.begin(), signal.end()));
242  }
243  if (wireReadoutGeom.SignalType(rdvec[rd]->Channel()) == geo::kCollection) {
244  if (*max_element(adc.begin(), adc.end()) > pedestal + threshold)
245  fRawColPeak->Fill(*max_element(adc.begin(), adc.end()));
246  if (*max_element(signal.begin(), signal.end()) > pedestal + threshold)
247  fCalColPeak->Fill(*max_element(signal.begin(), signal.end()));
248  }
249 
250  int window = 8;
251  static unsigned int pulseHeight = 5;
252  unsigned int tmin = 1;
253  unsigned int tmax = 1;
254  int indMax = TMath::LocMax(signalSize, &adc[0]);
255  double sigMin = 0.0;
256  double sigMax = TMath::MaxElement(signalSize, &adc[0]);
257  if (wireReadoutGeom.SignalType(rdvec[rd]->Channel()) == geo::kInduction &&
258  sigMax >= pulseHeight) {
259  int indMin = TMath::LocMin(signalSize, &adc[0]);
260  sigMin = TMath::MinElement(signalSize, &adc[0]);
261  tmin = std::max(indMax - window, 0);
262  tmax = std::min(indMin + window, (int)signalSize - 1);
263  MF_LOG_DEBUG("CalWireAna") << "Induction channel, indMin,tmin,tmax " << rd << " " << indMin
264  << " " << tmin << " " << tmax;
265  }
266  else if (sigMax >= pulseHeight) {
267  tmin = std::max(indMax - window, 0);
268  tmax = std::min(indMax + window, (int)signalSize - 1);
269  MF_LOG_DEBUG("CalWireAna")
270  << "Collection channel, tmin,tmax " << rd << " " << tmin << " " << tmax;
271  }
272 
273  fMin->Fill(sigMin);
274  fMax->Fill(sigMax);
275  fWindow->Fill(indMax, tmax - tmin);
276 
277  std::vector<double> winDiffs;
278  int cnt = 0;
279  static unsigned int tRawLead = 0;
280  for (unsigned int t = 1; t < signalSize; ++t) {
281  fDiffsW->Fill(signal[t] - signal[t - 1]);
282  fWireSig->Fill(rd, signal[t]);
283 
284  if (t >= tmin && t <= tmax && tmax >= pulseHeight && (t + tRawLead) < signalSize) {
285  cnt++;
286  winDiffs.push_back((adc[t + tRawLead] - signal[t]) / adc[t + tRawLead]);
287  fDiffsRW->Fill(adc[t + tRawLead] - signal[t]);
288  fDiffsRWvsR->Fill(rd, adc[t + tRawLead] - signal[t]);
289  if (sigMax >= pulseHeight) {
290  fDiffsRWgph->Fill(adc[t + tRawLead] - signal[t]);
291  fDiffsRWvsRgph->Fill(rd, adc[t + tRawLead] - signal[t]);
292  }
293  }
294  }
295 
296  MF_LOG_DEBUG("CalWireAna") << "on channel " << rdvec[rd]->Channel();
297  // TMath::Mean with iterators doesn't work. EC,23-Sep-2010.
298  double tmp = TMath::Mean(winDiffs.size(), &winDiffs[0]);
299  double tmp2 = TMath::RMS(winDiffs.size(), &winDiffs[0]);
300  for (int i = 0; i < fft->FFTSize(); i++) {
301  fIR->Fill(i, ir[i]);
302  fIW->Fill(i, iw[i]);
303  fCR->Fill(i, cr[i]);
304  fCW->Fill(i, cw[i]);
305  }
306  fRD_WireMeanDiff2D->Fill(rd, tmp);
307  fRD_WireRMSDiff2D->Fill(rd, tmp2);
308  } //end loop over rawDigits
309  } //end analyze method
310 
312 
313 } //end namespace
TH2F * fRD_WireMeanDiff2D
histogram of difference between original tdc value and compressesed value vs original value ...
STL namespace.
Definition of basic raw digits.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
void emplace_back(Args &&...args)
Definition: PtrVector.h:448
Float_t tmp
Definition: plot.C:35
void DoFFT(std::vector< T > &input, std::vector< TComplex > &output)
Definition: LArFFT.h:95
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
creation of calibrated signals on wires
TH1F * fDiffsR
histogram of Raw tdc to tdc differences
int FFTSize() const
Definition: LArFFT.h:66
TH2F * fRD_WireRMSDiff2D
histogram of difference between original tdc value and compressesed value vs original value ...
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
Signal from induction planes.
Definition: geo_types.h:147
void beginJob() override
size_type size() const
Definition: PtrVector.h:302
std::string fCalWireModuleLabel
name of module that produced the wires
CalWireAna(fhicl::ParameterSet const &pset)
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
void AlignedSum(std::vector< T > &input, std::vector< T > &output, bool add=true)
Definition: LArFFT.h:241
std::string fDetSimModuleLabel
#define MF_LOG_DEBUG(id)
TH1F * fDiffsW
histogram of Wire (post-deconvoution) tdc to tdc differences
Declaration of basic channel signal object.
TCEvent evt
Definition: DataStructs.cxx:8
Double_t sum
Definition: plot.C:31
void analyze(const art::Event &evt) override
read/write access to event
Base class for creation of raw signals on wires.
Signal from collection planes.
Definition: geo_types.h:148