148 = art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider();
149 art::Handle< std::vector<raw::RawDigit> > rdHandle;
151 art::Handle< std::vector<recob::Wire> > wHandle;
158 art::PtrVector<recob::Wire> wvec;
159 for(
unsigned int i = 0; i < wHandle->size(); ++i){
160 art::Ptr<recob::Wire>
w(wHandle,i);
163 art::PtrVector<raw::RawDigit> rdvec;
164 for(
unsigned int i = 0; i < rdHandle->size(); ++i){
165 art::Ptr<raw::RawDigit>
r(rdHandle,i);
169 art::ServiceHandle<geo::Geometry const> geom;
170 art::ServiceHandle<util::LArFFT> fft;
171 double pedestal = rdvec[0]->GetPedestal();
172 double threshold = 9.0;
173 double signalSize = rdvec[0]->Samples();
174 uint32_t indChan0=64;
175 uint32_t indChan1=110;
176 uint32_t colChan0=312;
177 uint32_t colChan1=354;
178 std::vector<double> ir(fft->FFTSize()),iw(fft->FFTSize()),cr(fft->FFTSize()),cw(fft->FFTSize());
179 ir[signalSize]=iw[signalSize]=cr[signalSize]=cw[signalSize]=1.0;
181 for(
unsigned int rd = 0; rd < rdvec.size(); ++rd){
183 std::vector<double> signal(fft->FFTSize());
184 for(
unsigned int wd = 0; wd < wvec.size(); ++wd){
186 if (RawDigitsFromWire.at(wd) == rdvec[rd]){
187 std::vector<float> wirSig = wvec[wd]->Signal();
188 if(wirSig.size() > signal.size()) {
189 MF_LOG_DEBUG(
"CalWireAna")<<
"Incompatible vector size "<<wirSig.size()
190 <<
" "<<signal.size();
193 for(
unsigned int ii = 0; ii < wirSig.size(); ++ii) {
194 signal[ii] = wirSig[ii];
198 if (wd == (wvec.size()-1) ){
199 MF_LOG_DEBUG(
"CalWireAna") <<
"caldata::CalWireAna:Big problem! No matching Wire for RawDigit. Bailing." << rd;
204 std::vector<double> adc(fft->FFTSize());
206 for(
unsigned int t = 1; t < rdvec[rd]->Samples(); ++t){
207 fDiffsR->Fill(rdvec[rd]->ADC(t) - rdvec[rd]->ADC(t-1));
208 adc[t-1]=rdvec[rd]->ADC(t-1);
209 fRawSig->Fill(rd,rdvec[rd]->ADC(t));
212 adc[rdvec[rd]->Samples()-1] = rdvec[rd]->ADC(rdvec[rd]->Samples()-1);
213 if(!channelStatus.
IsBad(rdvec[rd]->Channel()) &&
214 (*max_element(adc.begin(),adc.end()) < pedestal+threshold &&
215 *min_element(adc.begin(),adc.end()) >pedestal -threshold)) {
217 for(
int i = 0; i < signalSize; i++) sum+=pow(adc[i]-pedestal,2.0);
219 std::vector<double> temp(fft->FFTSize());
220 std::vector<TComplex> fTemp(fft->FFTSize()/2+1);
221 for(
int i = 0; i < signalSize; i++) temp[i]=(adc[i]-pedestal)*sin(TMath::Pi()*(
double)i/signalSize);
222 fft->DoFFT(temp,fTemp);
223 for(
int i = 0; i < fft->FFTSize()/2+1; i++)
fNoiseHist->Fill(i,fTemp[i].Rho());
226 rdvec[rd]->Channel() > indChan0 &&
227 rdvec[rd]->Channel() < indChan1){
228 fft->AlignedSum(ir,adc);
229 fft->AlignedSum(iw,signal);
232 rdvec[rd]->Channel() > colChan0 &&
233 rdvec[rd]->Channel() < colChan1) {
234 fft->AlignedSum(cr,adc);
235 fft->AlignedSum(cw,signal);
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()));
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()));
251 static unsigned int pulseHeight = 5;
252 unsigned int tmin = 1;
253 unsigned int tmax = 1;
254 int indMax = TMath::LocMax(signalSize,&adc[0]);
256 double sigMax = TMath::MaxElement(signalSize,&adc[0]);
257 if(geom->SignalType(rdvec[rd]->Channel()) ==
geo::kInduction && sigMax>=pulseHeight) {
258 int indMin = TMath::LocMin(signalSize,&adc[0]);
259 sigMin = TMath::MinElement(signalSize,&adc[0]);
260 tmin = std::max(indMax-window,0);
261 tmax = std::min(indMin+window,(
int)signalSize-1);
262 MF_LOG_DEBUG(
"CalWireAna") <<
"Induction channel, indMin,tmin,tmax "
263 << rd<<
" " << indMin<<
" " << tmin <<
" " << tmax;
265 else if (sigMax>=pulseHeight){
266 tmin = std::max(indMax-window,0);
267 tmax = std::min(indMax+window,(
int)signalSize-1);
268 MF_LOG_DEBUG(
"CalWireAna") <<
"Collection channel, tmin,tmax "<< rd<<
" " << tmin <<
" " << tmax;
273 fWindow->Fill(indMax, tmax - tmin);
275 std::vector<double> winDiffs;
278 static unsigned int tRawLead = 0;
279 for(
unsigned int t = 1; t < signalSize; ++t){
280 fDiffsW->Fill(signal[t]-signal[t-1]);
283 if (t>=tmin && t<=tmax && tmax>=pulseHeight && (t+tRawLead)<signalSize){
285 winDiffs.push_back((adc[t+tRawLead]-signal[t])/adc[t+tRawLead]);
286 fDiffsRW->Fill(adc[t+tRawLead]-signal[t]);
288 if (sigMax >= pulseHeight){
295 MF_LOG_DEBUG(
"CalWireAna") <<
"on channel " << rdvec[rd]->Channel();
297 double tmp = TMath::Mean(winDiffs.size(),&winDiffs[0]);
298 double tmp2 = TMath::RMS(winDiffs.size(),&winDiffs[0]);
300 for (
unsigned int ii=0; ii<rdvec[rd]->Samples(); ii++) tmp3+=rdvec[rd]->ADC(ii);
301 for(
int i = 0; i < fft->FFTSize(); i++) {
TH2F * fRD_WireMeanDiff2D
histogram of difference between original tdc value and compressesed value vs original value ...
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
TH1F * fDiffsR
histogram of Raw tdc to tdc differences
TH2F * fRD_WireRMSDiff2D
histogram of difference between original tdc value and compressesed value vs original value ...
Signal from induction planes.
Class providing information about the quality of channels.
std::string fCalWireModuleLabel
name of module that produced the wires
Var Sqrt(const Var &v)
Use to take sqrt of a var.
std::string fDetSimModuleLabel
TH1F * fDiffsW
histogram of Wire (post-deconvoution) tdc to tdc differences
Signal from collection planes.