Bug Summary

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