Bug Summary

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

Annotated Source Code

1// $Id$
2//
3// File: DPSHit_factory.cc
4// Created: Wed Oct 15 16:45:01 EDT 2014
5// Creator: staylor (on Linux gluon05.jlab.org 2.6.32-358.18.1.el6.x86_64 x86_64)
6//
7#include <iostream>
8#include <iomanip>
9using namespace std;
10
11#include "DPSHit_factory.h"
12#include <DAQ/Df250PulsePedestal.h>
13#include <DAQ/Df250PulseIntegral.h>
14#include <TTAB/DTTabUtilities.h>
15using namespace jana;
16
17
18//------------------
19// init
20//------------------
21jerror_t DPSHit_factory::init(void)
22{
23 ADC_THRESHOLD = 500.0; // ADC integral counts
24 gPARMS->SetDefaultParameter("PSHit:ADC_THRESHOLD",ADC_THRESHOLD,
25 "pedestal-subtracted pulse integral threshold");
26
27 /// set the base conversion scales
28 a_scale = 0.0001;
29 t_scale = 0.0625; // 62.5 ps/count
30 t_base = 0.; // ns
31
32 CHECK_FADC_ERRORS = false;
33 gPARMS->SetDefaultParameter("PSHit:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits");
34
35 return NOERROR;
36}
37
38//------------------
39// brun
40//------------------
41jerror_t DPSHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
42{
43 // Only print messages for one thread whenever run number change
44 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
45 static set<int> runs_announced;
46 pthread_mutex_lock(&print_mutex);
47 bool print_messages = false;
48 if(runs_announced.find(runnumber) == runs_announced.end()){
49 print_messages = true;
50 runs_announced.insert(runnumber);
51 }
52 pthread_mutex_unlock(&print_mutex);
53
54 /// Read in calibration constants
55 if(print_messages) jout << "In DPSHit_factory, loading constants..." << endl;
56
57 // extract the PS Geometry
58 vector<const DPSGeometry*> psGeomVect;
59 eventLoop->Get(psGeomVect);
60 if (psGeomVect.size() < 1)
61 return OBJECT_NOT_AVAILABLE;
62 const DPSGeometry& psGeom = *(psGeomVect[0]);
63
64 // load scale factors
65 map<string,double> scale_factors;
66 if (eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/digi_scales", scale_factors))
67 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/digi_scales !" << endl;
68 if (scale_factors.find("PS_ADC_ASCALE") != scale_factors.end())
69 a_scale = scale_factors["PS_ADC_ASCALE"];
70 else
71 jerr << "Unable to get PS_ADC_ASCALE from /PHOTON_BEAM/pair_spectrometer/digi_scales !"
72 << endl;
73 if (scale_factors.find("PS_ADC_TSCALE") != scale_factors.end())
74 t_scale = scale_factors["PS_ADC_TSCALE"];
75 else
76 jerr << "Unable to get PS_ADC_TSCALE from /PHOTON_BEAM/pair_spectrometer/digi_scales !"
77 << endl;
78
79 // load base time offset
80 map<string,double> base_time_offset;
81 if (eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/base_time_offset",base_time_offset))
82 jout << "Error loading /PHOTON_BEAM/pair_spectrometer/base_time_offset !" << endl;
83 if (base_time_offset.find("PS_FINE_BASE_TIME_OFFSET") != base_time_offset.end())
84 t_base = base_time_offset["PS_FINE_BASE_TIME_OFFSET"];
85 else
86 jerr << "Unable to get PS_FINE_BASE_TIME_OFFSET from /PHOTON_BEAM/pair_spectrometer/base_time_offset !" << endl;
87
88 FillCalibTable(adc_pedestals, "/PHOTON_BEAM/pair_spectrometer/fine/adc_pedestals", psGeom);
89 FillCalibTable(adc_gains, "/PHOTON_BEAM/pair_spectrometer/fine/adc_gain_factors", psGeom);
90 FillCalibTable(adc_time_offsets, "/PHOTON_BEAM/pair_spectrometer/fine/adc_timing_offsets", psGeom);
91
92 return NOERROR;
93}
94
95//------------------
96// evnt
97//------------------
98jerror_t DPSHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
99{
100 /// Generate DPSHit object for each DPSDigiHit object.
101 /// This is where the first set of calibration constants
102 /// is applied to convert from digitzed units into natural
103 /// units.
104 ///
105 /// Note that this code does NOT get called for simulated
106 /// data in HDDM format. The HDDM event source will copy
107 /// the precalibrated values directly into the _data vector.
108
109 // extract the PS Geometry
110 vector<const DPSGeometry*> psGeomVect;
111 eventLoop->Get(psGeomVect);
112 if (psGeomVect.size() < 1)
113 return OBJECT_NOT_AVAILABLE;
114 const DPSGeometry& psGeom = *(psGeomVect[0]);
115
116 const DTTabUtilities* locTTabUtilities = NULL__null;
117 loop->GetSingle(locTTabUtilities);
118
119 // First, make hits out of all fADC250 hits
120 vector<const DPSDigiHit*> digihits;
121 loop->Get(digihits);
122 char str[256];
123
124 for (unsigned int i=0; i < digihits.size(); i++){
125 const DPSDigiHit *digihit = digihits[i];
126
127 // Make sure channel id is in valid range
128 if( (digihit->arm < 0) || (digihit->arm >= psGeom.NUM_ARMS) ) {
129 sprintf(str, "DPSDigiHit arm out of range! arm=%d (should be 0-%d)",
130 digihit->arm, psGeom.NUM_ARMS);
131 throw JException(str);
132 }
133 if( (digihit->column <= 0) || (digihit->column > psGeom.NUM_FINE_COLUMNS) ) {
134 sprintf(str, "DPSDigiHit column out of range! column=%d (should be 0-%d)",
135 digihit->column, psGeom.NUM_FINE_COLUMNS);
136 throw JException(str);
137 }
138
139 // Throw away hits with firmware errors (post-summer 2016 firmware)
140 if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF))
141 continue;
142
143 if(digihit->datasource == 1) { // handle pre-Fall 2016 firmware
144 // Throw away hits where the fADC timing algorithm failed
145 //if (digihit->pulse_time == 0) continue;
146 // The following condition signals an error state in the flash algorithm
147 // Do not make hits out of these
148 const Df250PulsePedestal* PPobj = NULL__null;
149 digihit->GetSingle(PPobj);
150 if (PPobj != NULL__null){
151 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
152 }
153 else continue;
154 }
155
156 // Get pedestal, prefer associated event pedestal if it exists,
157 // otherwise, use the average pedestal from CCDB
158 double pedestal = GetConstant(adc_pedestals,digihit,psGeom);
Value stored to 'pedestal' during its initialization is never read
159 double nsamples_integral = (double)digihit->nsamples_integral;
160 double nsamples_pedestal = (double)digihit->nsamples_pedestal;
161
162 // nsamples_pedestal should always be positive for valid data - err on the side of caution for now
163 if(nsamples_pedestal == 0) {
164 jerr << "DPSDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl;
165 continue;
166 }
167
168 if ( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) {
169 // the measured pedestal must be scaled by the ratio of the number
170 // of samples used to calculate the integral and the pedestal
171 // Changed to conform to D. Lawrence changes Dec. 4 2014
172 double ped_sum = (double)digihit->pedestal;
173 pedestal = ped_sum * nsamples_integral/nsamples_pedestal;
174 }
175 else continue;
176
177 double single_sample_ped = pedestal/nsamples_pedestal;
178
179 // Apply calibration constants here
180 double P = digihit->pulse_peak - single_sample_ped;
181 double A = (double)digihit->pulse_integral;
182 A -= pedestal;
183 // Throw away hits with small pedestal-subtracted integrals
184 if (A < ADC_THRESHOLD) continue;
185 double T = (double)digihit->pulse_time;
186
187 DPSHit *hit = new DPSHit;
188 // The PSHit class labels hits as
189 // arm: North/South (0/1)
190 // column: 1-145
191 hit->arm = digihit->arm;
192 hit->column = digihit->column;
193 hit->integral = A;
194 hit->pulse_peak = P;
195 hit->npix_fadc = A * a_scale * GetConstant(adc_gains, digihit, psGeom);
196 hit->t = t_scale * T - GetConstant(adc_time_offsets, digihit, psGeom) + t_base;
197 hit->E = 0.5*(psGeom.getElow(digihit->arm,digihit->column) + psGeom.getEhigh(digihit->arm,digihit->column));
198
199 hit->AddAssociatedObject(digihit);
200
201 _data.push_back(hit);
202 }
203
204
205 return NOERROR;
206}
207
208//------------------
209// erun
210//------------------
211jerror_t DPSHit_factory::erun(void)
212{
213 return NOERROR;
214}
215
216//------------------
217// fini
218//------------------
219jerror_t DPSHit_factory::fini(void)
220{
221 return NOERROR;
222}
223
224//------------------
225// FillCalibTable
226//------------------
227void DPSHit_factory::FillCalibTable(ps_digi_constants_t &table, string table_name,
228 const DPSGeometry &psGeom)
229{
230 char str[256];
231
232 // load constant table
233 if(eventLoop->GetCalib(table_name, table))
234 jout << "Error loading " + table_name + " !" << endl;
235
236
237 // check that the size of the tables loaded are correct
238 if( (int)table.size() != psGeom.NUM_FINE_COLUMNS ) {
239 sprintf(str, "PS table loaded with wrong size! number of columns=%d (should be %d)",
240 (int)table.size(), psGeom.NUM_FINE_COLUMNS );
241 cerr << str << endl;
242 throw JException(str);
243 }
244
245 for( int column=0; column < psGeom.NUM_FINE_COLUMNS; column++) {
246 if( (int)table[column].size() != psGeom.NUM_ARMS ) {
247 sprintf(str, "PS table loaded with wrong size! column=%d number of arms=%d (should be %d)",
248 column, (int)table[column].size(), psGeom.NUM_ARMS );
249 cerr << str << endl;
250 throw JException(str);
251 }
252 }
253}
254
255//------------------------------------
256// GetConstant
257// Allow a few different interfaces
258//
259// PS Geometry as defined in the Translation Table:
260// arm: North/South (0/1)
261// column: 1-145
262// Note the different counting schemes used
263//------------------------------------
264const double DPSHit_factory::GetConstant( const ps_digi_constants_t &the_table,
265 const DPSGeometry::Arm in_arm, const int in_column,
266 const DPSGeometry &psGeom ) const
267{
268 char str[256];
269
270 if( (in_arm != DPSGeometry::kNorth) && (in_arm != DPSGeometry::kSouth)) {
271 sprintf(str, "Bad arm requested in DPSHit_factory::GetConstant()! requested=%d , should be 0-%d",
272 static_cast<int>(in_arm), static_cast<int>(DPSGeometry::kSouth));
273 cerr << str << endl;
274 throw JException(str);
275 }
276 if( (in_column <= 0) || (in_column > psGeom.NUM_FINE_COLUMNS)) {
277 sprintf(str, "Bad column # requested in DPSHit_factory::GetConstant()! requested=%d , should be 1-%d", in_column, psGeom.NUM_FINE_COLUMNS);
278 cerr << str << endl;
279 throw JException(str);
280 }
281
282 // the tables are indexed by column, with the different values for the two arms
283 // stored in the two fields of the pair
284 if(in_arm == DPSGeometry::kNorth) {
285 return the_table[in_column-1][in_arm];
286 } else {
287 return the_table[in_column-1][in_arm];
288 }
289}
290
291const double DPSHit_factory::GetConstant( const ps_digi_constants_t &the_table,
292 const DPSHit *in_hit, const DPSGeometry &psGeom ) const
293{
294 char str[256];
295
296 if( (in_hit->arm != DPSGeometry::kNorth) && (in_hit->arm != DPSGeometry::kSouth)) {
297 sprintf(str, "Bad arm requested in DPSHit_factory::GetConstant()! requested=%d , should be 0-%d",
298 static_cast<int>(in_hit->arm), static_cast<int>(DPSGeometry::kSouth));
299 cerr << str << endl;
300 throw JException(str);
301 }
302 if( (in_hit->column <= 0) || (in_hit->column > psGeom.NUM_FINE_COLUMNS)) {
303 sprintf(str, "Bad column # requested in DPSHit_factory::GetConstant()! requested=%d , should be 1-%d", in_hit->column, psGeom.NUM_FINE_COLUMNS);
304 cerr << str << endl;
305 throw JException(str);
306 }
307
308 // the tables are indexed by column, with the different values for the two arms
309 // stored in the two fields of the pair
310 if(in_hit->arm == DPSGeometry::kNorth) {
311 return the_table[in_hit->column-1][in_hit->arm];
312 } else {
313 return the_table[in_hit->column-1][in_hit->arm];
314 }
315}
316
317const double DPSHit_factory::GetConstant( const ps_digi_constants_t &the_table,
318 const DPSDigiHit *in_digihit, const DPSGeometry &psGeom) const
319{
320 char str[256];
321
322 if( (in_digihit->arm != DPSGeometry::kNorth) && (in_digihit->arm != DPSGeometry::kSouth)) {
323 sprintf(str, "Bad arm requested in DPSHit_factory::GetConstant()! requested=%d , should be 0-%d",
324 static_cast<int>(in_digihit->arm), static_cast<int>(DPSGeometry::kSouth));
325 cerr << str << endl;
326 throw JException(str);
327 }
328 if( (in_digihit->column <= 0) || (in_digihit->column > psGeom.NUM_FINE_COLUMNS)) {
329 sprintf(str, "Bad column # requested in DPSHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->column, psGeom.NUM_FINE_COLUMNS);
330 cerr << str << endl;
331 throw JException(str);
332 }
333
334 // the tables are indexed by column, with the different values for the two arms
335 // stored in the two fields of the pair
336 if(in_digihit->arm == DPSGeometry::kNorth) {
337 return the_table[in_digihit->column-1][in_digihit->arm];
338 } else {
339 return the_table[in_digihit->column-1][in_digihit->arm];
340 }
341}
342