Bug Summary

File:libraries/PAIR_SPECTROMETER/DPSCHit_factory.cc
Location:line 191, column 12
Description:Value stored to 'pedestal' during its initialization is never read

Annotated Source Code

1// $Id$
2//
3// File: DPSCHit_factory.cc
4// Created: Wed Oct 15 16:45:33 EDT 2014
5// Creator: staylor (on Linux gluon05.jlab.org 2.6.32-358.18.1.el6.x86_64 x86_64)
6//
7
8#include <iostream>
9#include <iomanip>
10#include <cmath>
11#include <vector>
12#include <limits>
13using namespace std;
14
15#include "DPSCHit_factory.h"
16#include "DPSCDigiHit.h"
17#include "DPSCTDCDigiHit.h"
18#include <DAQ/Df250PulsePedestal.h>
19#include <DAQ/Df250PulseIntegral.h>
20using namespace jana;
21
22
23//------------------
24// init
25//------------------
26jerror_t DPSCHit_factory::init(void)
27{
28 DELTA_T_ADC_TDC_MAX = 4.0; // ns
29 gPARMS->SetDefaultParameter("PSCHit:DELTA_T_ADC_TDC_MAX", DELTA_T_ADC_TDC_MAX,
30 "Maximum difference in ns between a (calibrated) fADC time and"
31 " F1TDC time for them to be matched in a single hit");
32 ADC_THRESHOLD = 500.0; // ADC integral counts
33 gPARMS->SetDefaultParameter("PSCHit:ADC_THRESHOLD",ADC_THRESHOLD,
34 "pedestal-subtracted pulse integral threshold");
35
36 CHECK_FADC_ERRORS = false;
37 gPARMS->SetDefaultParameter("PSCHit:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits");
38
39 /// set the base conversion scales
40 a_scale = 0.0001;
41 t_scale = 0.0625; // 62.5 ps/count
42 t_base = 0.; // ns
43 t_tdc_base = 0.;
44
45 return NOERROR;
46}
47
48//------------------
49// brun
50//------------------
51jerror_t DPSCHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
52{
53 // Only print messages for one thread whenever run number change
54 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
55 static set<int> runs_announced;
56 pthread_mutex_lock(&print_mutex);
57 bool print_messages = false;
58 if(runs_announced.find(runnumber) == runs_announced.end()){
59 print_messages = true;
60 runs_announced.insert(runnumber);
61 }
62 pthread_mutex_unlock(&print_mutex);
63
64 /// Read in calibration constants
65 if(print_messages) jout << "In DPSCHit_factory, loading constants..." << endl;
66
67 // extract the PS Geometry
68 vector<const DPSGeometry*> psGeomVect;
69 eventLoop->Get( psGeomVect );
70 if (psGeomVect.size() < 1)
71 return OBJECT_NOT_AVAILABLE;
72 const DPSGeometry& psGeom = *(psGeomVect[0]);
73
74 // load scale factors
75 map<string,double> scale_factors;
76 if (eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/digi_scales", scale_factors))
77 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/digi_scales !" << endl;
78 if (scale_factors.find("PSC_ADC_ASCALE") != scale_factors.end())
79 a_scale = scale_factors["PSC_ADC_ASCALE"];
80 else
81 jerr << "Unable to get PSC_ADC_ASCALE from /PHOTON_BEAM/pair_spectrometer/digi_scales !"
82 << endl;
83 if (scale_factors.find("PSC_ADC_TSCALE") != scale_factors.end())
84 t_scale = scale_factors["PSC_ADC_TSCALE"];
85 else
86 jerr << "Unable to get PSC_ADC_TSCALE from /PHOTON_BEAM/pair_spectrometer/digi_scales !"
87 << endl;
88
89 // load base time offset
90 map<string,double> base_time_offset;
91 if (eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/base_time_offset",base_time_offset))
92 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/base_time_offset !" << endl;
93 if (base_time_offset.find("PS_COARSE_BASE_TIME_OFFSET") != base_time_offset.end())
94 t_base = base_time_offset["PS_COARSE_BASE_TIME_OFFSET"];
95 else
96 jerr << "Unable to get PS_COARSE_BASE_TIME_OFFSET from /PHOTON_BEAM/pair_spectrometer/base_time_offset !" << endl;
97 //
98 if (base_time_offset.find("PS_COARSE_TDC_BASE_TIME_OFFSET") != base_time_offset.end())
99 t_tdc_base = base_time_offset["PS_COARSE_TDC_BASE_TIME_OFFSET"];
100 else
101 jerr << "Unable to get PS_COARSE_TDC_BASE_TIME_OFFSET from /PHOTON_BEAM/pair_spectrometer/base_time_offset !" << endl;
102
103 /// Read in calibration constants
104 vector<double> raw_adc_pedestals;
105 vector<double> raw_adc_gains;
106 vector<double> raw_adc_offsets;
107 vector<double> raw_tdc_offsets;
108
109 // load constant tables
110 if(eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/coarse/adc_pedestals", raw_adc_pedestals))
111 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/coarse/adc_pedestals !" << endl;
112 if(eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/coarse/adc_gain_factors", raw_adc_gains))
113 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/coarse/adc_gain_factors !" << endl;
114 if(eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/coarse/adc_timing_offsets", raw_adc_offsets))
115 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/coarse/adc_timing_offsets !" << endl;
116 if(eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/coarse/tdc_timing_offsets", raw_tdc_offsets))
117 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/coarse/tdc_timing_offsets !" << endl;
118
119
120 FillCalibTable(adc_pedestals, raw_adc_pedestals, psGeom);
121 FillCalibTable(adc_gains, raw_adc_gains, psGeom);
122 FillCalibTable(adc_time_offsets, raw_adc_offsets, psGeom);
123 FillCalibTable(tdc_time_offsets, raw_tdc_offsets, psGeom);
124
125 // load timewalk corrections
126 if (eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/tdc_timewalk_corrections", tw_parameters))
127 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/tdc_timewalk_corrections !" << endl;
128
129 return NOERROR;
130}
131
132//------------------
133// evnt
134//------------------
135jerror_t DPSCHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
136{
137 /// Generate DPSCHit object for each DPSCDigiHit object.
138 /// This is where the first set of calibration constants
139 /// is applied to convert from digitzed units into natural
140 /// units.
141 ///
142 /// Note that this code does NOT get called for simulated
143 /// data in HDDM format. The HDDM event source will copy
144 /// the precalibrated values directly into the _data vector.
145
146 // extract the PS Geometry
147 vector<const DPSGeometry*> psGeomVect;
148 eventLoop->Get(psGeomVect);
149 if (psGeomVect.size() < 1)
150 return OBJECT_NOT_AVAILABLE;
151 const DPSGeometry& psGeom = *(psGeomVect[0]);
152
153 const DTTabUtilities* locTTabUtilities = NULL__null;
154 loop->GetSingle(locTTabUtilities);
155
156 // First, make hits out of all fADC250 hits
157 vector<const DPSCDigiHit*> digihits;
158 loop->Get(digihits);
159 char str[256];
160
161 for (unsigned int i=0; i < digihits.size(); i++){
162 const DPSCDigiHit *digihit = digihits[i];
163
164 // Make sure channel id is in valid range
165 const int PSC_MAX_CHANNELS = psGeom.NUM_COARSE_COLUMNS*psGeom.NUM_ARMS;
166 if( (digihit->counter_id <= 0) || (digihit->counter_id > PSC_MAX_CHANNELS)) {
167 sprintf(str, "DPSCDigiHit sector out of range! sector=%d (should be 1-%d)",
168 digihit->counter_id, PSC_MAX_CHANNELS);
169 throw JException(str);
170 }
171
172 // Throw away hits with firmware errors (post-summer 2016 firmware)
173 if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF))
174 continue;
175
176 if(digihit->datasource == 1) { // handle pre-Fall 2016 firmware
177 // Throw away hits where the fADC timing algorithm failed
178 //if (digihit->pulse_time == 0) continue;
179 // The following condition signals an error state in the flash algorithm
180 // Do not make hits out of these
181 const Df250PulsePedestal* PPobj = NULL__null;
182 digihit->GetSingle(PPobj);
183 if (PPobj != NULL__null){
184 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
185 }
186 else continue;
187 }
188
189 // Get pedestal, prefer associated event pedestal if it exists,
190 // otherwise, use the average pedestal from CCDB
191 double pedestal = GetConstant(adc_pedestals,digihit,psGeom);
Value stored to 'pedestal' during its initialization is never read
192 double nsamples_integral = (double)digihit->nsamples_integral;
193 double nsamples_pedestal = (double)digihit->nsamples_pedestal;
194
195 // nsamples_pedestal should always be positive for valid data - err on the side of caution for now
196 if(nsamples_pedestal == 0) {
197 jerr << "DPSCDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl;
198 continue;
199 }
200
201 if ( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) {
202 // the measured pedestal must be scaled by the ratio of the number
203 // of samples used to calculate the integral and the pedestal
204 // Changed to conform to D. Lawrence changes Dec. 4 2014
205 double ped_sum = (double)digihit->pedestal;
206 pedestal = ped_sum * nsamples_integral/nsamples_pedestal;
207 }
208 else continue;
209
210 double single_sample_ped = pedestal/nsamples_pedestal;
211
212 // Apply calibration constants here
213 double P = digihit->pulse_peak - single_sample_ped;
214 double A = (double)digihit->pulse_integral;
215 A -= pedestal;
216 // Throw away hits with small pedestal-subtracted integrals
217 if (A < ADC_THRESHOLD) continue;
218 double T = (double)digihit->pulse_time;
219
220 DPSCHit *hit = new DPSCHit;
221 hit->arm = GetArm(digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
222 hit->module = GetModule(digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
223 hit->integral = A;
224 hit->pulse_peak = P;
225 hit->npe_fadc = A * a_scale * GetConstant(adc_gains, digihit, psGeom);
226 hit->time_fadc = t_scale * T - GetConstant(adc_time_offsets, digihit, psGeom) + t_base;
227 hit->t = hit->time_fadc;
228 hit->time_tdc = numeric_limits<double>::quiet_NaN();
229 hit->has_fADC = true;
230 hit->has_TDC = false; // will get set to true below if appropriate
231
232 hit->AddAssociatedObject(digihit);
233
234 _data.push_back(hit);
235 }
236
237 // Second, loop over TDC hits, matching them to the
238 // existing fADC hits where possible and updating
239 // their time information. If no match is found, then
240 // create a new hit with just the TDC info.
241 vector<const DPSCTDCDigiHit*> tdcdigihits;
242 loop->Get(tdcdigihits);
243
244 for(unsigned int i=0; i<tdcdigihits.size(); i++) {
245 const DPSCTDCDigiHit *digihit = tdcdigihits[i];
246
247 // calculate geometry information
248 DPSGeometry::Arm arm = GetArm(digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
249 int module = GetModule(digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
250
251 // Apply calibration constants here
252 double T = locTTabUtilities->Convert_DigiTimeToNs_F1TDC(digihit) - GetConstant(tdc_time_offsets, digihit, psGeom) + t_tdc_base;
253 // Look for existing hits to see if there is a match
254 // or create new one if there is no match
255 DPSCHit *hit = FindMatch(arm, module, T);
256 if (!hit) {
257 hit = new DPSCHit;
258 hit->arm = arm;
259 hit->module = module;
260 hit->time_fadc = numeric_limits<double>::quiet_NaN();
261 hit->integral = numeric_limits<double>::quiet_NaN();
262 hit->pulse_peak = numeric_limits<double>::quiet_NaN();
263 hit->npe_fadc = numeric_limits<double>::quiet_NaN();
264 hit->has_fADC = false;
265 _data.push_back(hit);
266 }
267 hit->time_tdc = T;
268 hit->has_TDC = true;
269 // apply time-walk corrections
270 if (hit->pulse_peak > 0) {
271 if (tw_parameters[module - 1][0] == 0) {
272 c1 = tw_parameters[module - 1][3];
273 c2 = tw_parameters[module - 1][4];
274 thresh = tw_parameters[module - 1][5];
275 P_0 = tw_parameters[module - 1][6];
276 }
277 else {
278 c1 = tw_parameters[module + 7][3];
279 c2 = tw_parameters[module + 7][4];
280 thresh = tw_parameters[module + 7][5];
281 P_0 = tw_parameters[module + 7][6];
282 }
283 double P = hit->pulse_peak;
284 T -= c1*(pow(P/thresh,c2) - pow(P_0/thresh,c2));
285 }
286 hit->t = T;
287
288 hit->AddAssociatedObject(digihit);
289 }
290
291 return NOERROR;
292}
293
294//------------------
295// FindMatch
296//------------------
297DPSCHit* DPSCHit_factory::FindMatch(DPSGeometry::Arm arm, int module, double T)
298{
299 DPSCHit* best_match = NULL__null;
300
301 // Loop over existing hits (from fADC) and look for a match
302 // in both the sector and the time.
303 for(unsigned int i=0; i<_data.size(); i++) {
304 DPSCHit *hit = _data[i];
305
306 if(!hit->has_fADC) continue; // only match to fADC hits, not bachelor TDC hits
307 if(hit->arm != arm) continue;
308 if(hit->module != module) continue;
309
310 double delta_T = fabs(T - hit->t);
311 if(delta_T > DELTA_T_ADC_TDC_MAX) continue;
312
313 return hit;
314
315 // if there are multiple hits, pick the one that is closest in time
316 if(best_match != NULL__null) {
317 if(delta_T < fabs(best_match->t - T))
318 best_match = hit;
319 } else {
320 best_match = hit;
321 }
322
323 }
324
325 return best_match;
326}
327
328// The translation table has PSC channels labaled as paddles 1-16
329// The PSCHit class labels hits as
330// arm: North/South (0/1)
331// module: 1-8
332const DPSGeometry::Arm DPSCHit_factory::GetArm(const int counter_id,const int num_counters_per_arm) const {
333 return counter_id <= num_counters_per_arm ? DPSGeometry::kNorth : DPSGeometry::kSouth;
334}
335const int DPSCHit_factory::GetModule(const int counter_id,const int num_counters_per_arm) const {
336 return counter_id <= num_counters_per_arm ? counter_id : counter_id - num_counters_per_arm;
337}
338
339//------------------
340// erun
341//------------------
342jerror_t DPSCHit_factory::erun(void)
343{
344 return NOERROR;
345}
346
347//------------------
348// fini
349//------------------
350jerror_t DPSCHit_factory::fini(void)
351{
352 return NOERROR;
353}
354
355//------------------
356// FillCalibTable
357//------------------
358void DPSCHit_factory::FillCalibTable(psc_digi_constants_t &table, vector<double> &raw_table,
359 const DPSGeometry &psGeom)
360{
361 char str[256];
362 const int PSC_MAX_CHANNELS = psGeom.NUM_COARSE_COLUMNS*psGeom.NUM_ARMS;
363 int channel = 0;
364
365 // reset the table before filling it
366 table.clear();
367
368 // initialize table
369 for(int column=0; column<psGeom.NUM_COARSE_COLUMNS; column++)
370 table.push_back( pair<double,double>() );
371
372 // the constants are stored in the CCDB as a 1D array
373 // with the first 8 channels being the north arm,
374 // and the second 8 channels being the south arm
375 for(int column=0; column<psGeom.NUM_COARSE_COLUMNS; column++) {
376 if( column+psGeom.NUM_COARSE_COLUMNS > PSC_MAX_CHANNELS) { // sanity check
377 sprintf(str, "Too many channels for PSC table! channel=%d (should be %d)",
378 channel, PSC_MAX_CHANNELS);
379 cerr << str << endl;
380 throw JException(str);
381 }
382
383 table[column] = pair<double,double>(raw_table[column],raw_table[column+psGeom.NUM_COARSE_COLUMNS]);
384 channel += 2;
385 }
386
387 // check to make sure that we loaded enough channels
388 if(channel != PSC_MAX_CHANNELS) {
389 sprintf(str, "Not enough channels for PSC table! channel=%d (should be %d)",
390 channel, PSC_MAX_CHANNELS);
391 cerr << str << endl;
392 throw JException(str);
393 }
394}
395
396//------------------------------------
397// GetConstant
398// Allow a few different interfaces
399//
400// PSC Geometry as defined in the Translation Table:
401// arm: North/South (0/1)
402// module: 1-8
403// Note the different counting schemes used
404//------------------------------------
405const double DPSCHit_factory::GetConstant( const psc_digi_constants_t &the_table,
406 const DPSGeometry::Arm in_arm, const int in_module,
407 const DPSGeometry &psGeom ) const
408{
409 char str[256];
410
411 if( (in_arm != DPSGeometry::kNorth) && (in_arm != DPSGeometry::kSouth)) {
412 sprintf(str, "Bad arm requested in DPSCHit_factory::GetConstant()! requested=%d , should be 0-%d",
413 static_cast<int>(in_arm), static_cast<int>(DPSGeometry::kSouth));
414 cerr << str << endl;
415 throw JException(str);
416 }
417 if( (in_module <= 0) || (in_module > psGeom.NUM_COARSE_COLUMNS)) {
418 sprintf(str, "Bad module # requested in DPSCHit_factory::GetConstant()! requested=%d , should be 1-%d", in_module, psGeom.NUM_COARSE_COLUMNS);
419 cerr << str << endl;
420 throw JException(str);
421 }
422
423 // the tables are indexed by module, with the different values for the two arms
424 // stored in the two fields of the pair
425 if(in_arm == DPSGeometry::kNorth) {
426 return the_table[in_module-1].first;
427 } else {
428 return the_table[in_module-1].second;
429 }
430}
431
432const double DPSCHit_factory::GetConstant( const psc_digi_constants_t &the_table,
433 const DPSCHit *in_hit, const DPSGeometry &psGeom ) const
434{
435 char str[256];
436
437 if( (in_hit->arm != DPSGeometry::kNorth) && (in_hit->arm != DPSGeometry::kSouth)) {
438 sprintf(str, "Bad arm requested in DPSCHit_factory::GetConstant()! requested=%d , should be 0-%d",
439 static_cast<int>(in_hit->arm), static_cast<int>(DPSGeometry::kSouth));
440 cerr << str << endl;
441 throw JException(str);
442 }
443 if( (in_hit->module <= 0) || (in_hit->module > psGeom.NUM_COARSE_COLUMNS)) {
444 sprintf(str, "Bad module # requested in DPSCHit_factory::GetConstant()! requested=%d , should be 1-%d", in_hit->module, psGeom.NUM_COARSE_COLUMNS);
445 cerr << str << endl;
446 throw JException(str);
447 }
448
449 // the tables are indexed by module, with the different values for the two arms
450 // stored in the two fields of the pair
451 if(in_hit->arm == DPSGeometry::kNorth) {
452 return the_table[in_hit->module-1].first;
453 } else {
454 return the_table[in_hit->module-1].second;
455 }
456}
457
458const double DPSCHit_factory::GetConstant( const psc_digi_constants_t &the_table,
459 const DPSCDigiHit *in_digihit, const DPSGeometry &psGeom) const
460{
461 char str[256];
462 // calculate geometry information
463 DPSGeometry::Arm arm = GetArm(in_digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
464 int module = GetModule(in_digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
465
466 if( (module <= 0) || (module > psGeom.NUM_COARSE_COLUMNS)) {
467 sprintf(str, "Bad module # requested in DPSCHit_factory::GetConstant()! requested=%d , should be 1-%d", module, psGeom.NUM_COARSE_COLUMNS);
468 cerr << str << endl;
469 throw JException(str);
470 }
471
472 // the tables are indexed by module, with the different values for the two arms
473 // stored in the two fields of the pair
474 if(arm == DPSGeometry::kNorth) {
475 return the_table[module-1].first;
476 } else {
477 return the_table[module-1].second;
478 }
479}
480
481const double DPSCHit_factory::GetConstant( const psc_digi_constants_t &the_table,
482 const DPSCTDCDigiHit *in_digihit, const DPSGeometry &psGeom ) const
483{
484 char str[256];
485 // calculate geometry information
486 DPSGeometry::Arm arm = GetArm(in_digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
487 int module = GetModule(in_digihit->counter_id,psGeom.NUM_COARSE_COLUMNS);
488
489 if( (module <= 0) || (module > psGeom.NUM_COARSE_COLUMNS)) {
490 sprintf(str, "Bad module # requested in DPSCHit_factory::GetConstant()! requested=%d , should be 1-%d", module, psGeom.NUM_COARSE_COLUMNS);
491 cerr << str << endl;
492 throw JException(str);
493 }
494
495 // the tables are indexed by module, with the different values for the two arms
496 // stored in the two fields of the pair
497 if(arm == DPSGeometry::kNorth) {
498 return the_table[module-1].first;
499 } else {
500 return the_table[module-1].second;
501 }
502}
503