Bug Summary

File:libraries/START_COUNTER/DSCHit_factory.cc
Location:line 201, column 16
Description:Value stored to 'pedestal' during its initialization is never read

Annotated Source Code

1// $Id$
2//
3// File: DSCHit_factory.cc
4// Created: Tue Aug 6 12:53:32 EDT 2013
5// Creator: davidl (on Darwin harriet.jlab.org 11.4.2 i386)
6//
7
8
9#include <iostream>
10#include <iomanip>
11#include <limits>
12#include <cmath>
13using namespace std;
14
15#include <START_COUNTER/DSCDigiHit.h>
16#include <START_COUNTER/DSCTDCDigiHit.h>
17#include <DAQ/Df250PulsePedestal.h>
18#include <DAQ/Df250PulseIntegral.h>
19#include <DAQ/Df250Config.h>
20#include "DSCHit_factory.h"
21
22using namespace jana;
23
24bool DSCHit_fadc_cmp(const DSCDigiHit *a,const DSCDigiHit *b)
25{
26 if (a->sector==b->sector) return (a->pulse_time<b->pulse_time);
27 return (a->sector<b->sector);
28}
29
30bool DSCHit_tdc_cmp(const DSCTDCDigiHit *a,const DSCTDCDigiHit *b)
31{
32 if (a->sector==b->sector) return (a->time<b->time);
33 return (a->sector<b->sector);
34}
35
36//------------------
37// init
38//------------------
39jerror_t DSCHit_factory::init(void)
40{
41 DELTA_T_ADC_TDC_MAX = 20.0; // ns
42 //DELTA_T_ADC_TDC_MAX = 50.0; // ns
43 //DELTA_T_ADC_TDC_MAX = 3600.0; // ns
44 gPARMS->SetDefaultParameter("SC:DELTA_T_ADC_TDC_MAX", DELTA_T_ADC_TDC_MAX,
45 "Maximum difference in ns between a (calibrated) fADC time and"
46 " F1TDC time for them to be matched in a single hit");
47
48 HIT_TIME_WINDOW = 60.0; //ns
49 gPARMS->SetDefaultParameter("SC:HIT_TIME_WINDOW", HIT_TIME_WINDOW,
50 "Time window of trigger corrected TDC time in which a hit in"
51 " in the TDC will match to a hit in the fADC to form an ST hit");
52
53 //ADC_THRESHOLD = 200.; // adc counts (= 50 mV threshold)
54 ADC_THRESHOLD = 120.; // adc counts (= 10 Mv threshold)
55 gPARMS->SetDefaultParameter("SC:ADC_THRESHOLD", ADC_THRESHOLD,
56 "Software pulse integral threshold");
57
58 USE_TIMEWALK_CORRECTION = 1.;
59 gPARMS->SetDefaultParameter("SC:USE_TIMEWALK_CORRECTION", USE_TIMEWALK_CORRECTION,
60 "Flag to decide if timewalk corrections should be applied.");
61
62 /// set the base conversion scales
63 a_scale = 0.0001;
64 t_scale = 0.0625; // 62.5 ps/count
65 t_base = 0.; // ns
66 t_tdc_base = 0.;
67
68 CHECK_FADC_ERRORS = false;
69 gPARMS->SetDefaultParameter("SC:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits");
70
71 return NOERROR;
72}
73
74//------------------
75// brun
76//------------------
77jerror_t DSCHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
78{
79 // Only print messages for one thread whenever run number change
80 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
81 static set<int> runs_announced;
82 pthread_mutex_lock(&print_mutex);
83 bool print_messages = false;
84 if(runs_announced.find(runnumber) == runs_announced.end()){
85 print_messages = true;
86 runs_announced.insert(runnumber);
87 }
88 pthread_mutex_unlock(&print_mutex);
89
90 /// Read in calibration constants
91 if(print_messages) jout << "In DSCHit_factory, loading constants..." << endl;
92
93 // load scale factors
94 map<string,double> scale_factors;
95 // a_scale (SC_ADC_SCALE)
96 if (eventLoop->GetCalib("/START_COUNTER/digi_scales", scale_factors))
97 jout << "Error loading /START_COUNTER/digi_scales !" << endl;
98 if (scale_factors.find("SC_ADC_ASCALE") != scale_factors.end())
99 a_scale = scale_factors["SC_ADC_ASCALE"];
100 else
101 jerr << "Unable to get SC_ADC_ASCALE from /START_COUNTER/digi_scales !"
102 << endl;
103 // t_scale (SC_ADC_SCALE)
104 if (scale_factors.find("SC_ADC_TSCALE") != scale_factors.end())
105 t_scale = scale_factors["SC_ADC_TSCALE"];
106 else
107 jerr << "Unable to get SC_ADC_TSCALE from /START_COUNTER/digi_scales !"
108 << endl;
109
110 // load base time offset
111 map<string,double> base_time_offset;
112 // t_base (SC_BASE_TIME_OFFSET)
113 if (eventLoop->GetCalib("/START_COUNTER/base_time_offset",base_time_offset))
114 jout << "Error loading /START_COUNTER/base_time_offset !" << endl;
115 if (base_time_offset.find("SC_BASE_TIME_OFFSET") != base_time_offset.end())
116 t_base = base_time_offset["SC_BASE_TIME_OFFSET"];
117 else
118 jerr << "Unable to get SC_BASE_TIME_OFFSET from /START_COUNTER/base_time_offset !" << endl;
119 // t_tdc_base (SC_TDC_BASE_TIME_OFFSET)
120 if (base_time_offset.find("SC_TDC_BASE_TIME_OFFSET") != base_time_offset.end())
121 t_tdc_base = base_time_offset["SC_TDC_BASE_TIME_OFFSET"];
122 else
123 jerr << "Unable to get SC_BASE_TIME_OFFSET from /START_COUNTER/base_time_offset !" << endl;
124
125 // load constant tables
126 // a_gains (gains)
127 if (eventLoop->GetCalib("/START_COUNTER/gains", a_gains))
128 jout << "Error loading /START_COUNTER/gains !" << endl;
129 // a_pedestals (pedestals)
130 if (eventLoop->GetCalib("/START_COUNTER/pedestals", a_pedestals))
131 jout << "Error loading /START_COUNTER/pedestals !" << endl;
132 // adc_time_offsets (adc_timing_offsets)
133 if (eventLoop->GetCalib("/START_COUNTER/adc_timing_offsets", adc_time_offsets))
134 jout << "Error loading /START_COUNTER/adc_timing_offsets !" << endl;
135 // tdc_time_offsets (tdc_timing_offsets)
136 if (eventLoop->GetCalib("/START_COUNTER/tdc_timing_offsets", tdc_time_offsets))
137 jout << "Error loading /START_COUNTER/tdc_timing_offsets !" << endl;
138 // timewalk_parameters (timewalk_parms)
139 if(eventLoop->GetCalib("START_COUNTER/timewalk_parms_v2", timewalk_parameters))
140 jout << "Error loading /START_COUNTER/timewalk_parms_v2 !" << endl;
141
142
143 return NOERROR;
144}
145
146//------------------
147// evnt
148//------------------
149jerror_t DSCHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
150{
151 /// Generate DSCHit object for each DSCDigiHit object.
152 /// This is where the first set of calibration constants
153 /// is applied to convert from digitzed units into natural
154 /// units.
155 ///
156 /// Note that this code does NOT get called for simulated
157 /// data in HDDM format. The HDDM event source will copy
158 /// the precalibrated values directly into the _data vector.
159
160 // First, make hits out of all fADC250 hits
161 vector<const DSCDigiHit*> digihits;
162 loop->Get(digihits);
163 sort(digihits.begin(),digihits.end(),DSCHit_fadc_cmp);
164
165 const DTTabUtilities* locTTabUtilities = NULL__null;
166 loop->GetSingle(locTTabUtilities);
167
168 char str[256];
169
170 for (unsigned int i = 0; i < digihits.size(); i++) {
171 const DSCDigiHit *digihit = digihits[i];
172
173 // Make sure sector is in valid range
174 if( (digihit->sector <= 0) && (digihit->sector > MAX_SECTORS)) {
175 sprintf(str, "DSCDigiHit sector out of range! sector=%d (should be 1-%d)",
176 digihit->sector, MAX_SECTORS);
177 throw JException(str);
178 }
179
180 // Error checking for pre-Fall 2016 firmware
181 if(digihit->datasource == 1) {
182 // There is a slight difference between Mode 7 and 8 data
183 // The following condition signals an error state in the flash algorithm
184 // Do not make hits out of these
185 const Df250PulsePedestal* PPobj = NULL__null;
186 digihit->GetSingle(PPobj);
187 if (PPobj != NULL__null) {
188 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
189 } else
190 continue;
191
192 if (digihit->pulse_time == 0) continue; // Should already be caught, but I'll leave it
193 }
194
195 if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF))
196 continue;
197
198 // Initialize pedestal to one found in CCDB, but override it
199 // with one found in event if is available (?)
200 // For now, only keep events with a correct pedestal
201 double pedestal = a_pedestals[digihit->sector-1];
Value stored to 'pedestal' during its initialization is never read
202 double nsamples_integral = digihit->nsamples_integral;
203 double nsamples_pedestal = digihit->nsamples_pedestal;
204
205 // nsamples_pedestal should always be positive for valid data - err on the side of caution for now
206 if(nsamples_pedestal == 0) {
207 jerr << "DSCDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl;
208 continue;
209 }
210
211 if( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) {
212 pedestal = digihit->pedestal * (nsamples_integral/nsamples_pedestal);
213 } else {
214 continue;
215 }
216
217 double single_sample_ped = pedestal/nsamples_pedestal;
218
219 //if ( ((double)digihit->pulse_integral) < ADC_THRESHOLD) continue; // Will comment out until this is set to something useful by default
220
221 double A = (double)digihit->pulse_integral;
222 double T = (double)digihit->pulse_time;
223
224 DSCHit *hit = new DSCHit;
225 // Sectors are numbered from 1-30
226 hit->sector = digihit->sector;
227
228 hit->dE = a_scale * a_gains[digihit->sector-1] * (A - pedestal);
229 hit->t_fADC = t_scale * T - adc_time_offsets[hit->sector-1] + t_base;
230 hit->t_TDC = numeric_limits<double>::quiet_NaN();
231
232 hit->has_TDC = false;
233 hit->has_fADC = true;
234
235 hit->t = hit->t_fADC; // set time from fADC in case no TDC hit
236 hit->pulse_height = digihit->pulse_peak - single_sample_ped;
237
238 hit->AddAssociatedObject(digihit);
239
240 _data.push_back(hit);
241 }
242
243
244 // Next, loop over TDC hits, matching them to the
245 // existing fADC hits where possible and updating
246 // their time information. If no match is found, then
247 // create a new hit with just the TDC info.
248 vector<const DSCTDCDigiHit*> tdcdigihits;
249 loop->Get(tdcdigihits);
250 sort(tdcdigihits.begin(),tdcdigihits.end(),DSCHit_tdc_cmp);
251
252 for (unsigned int i = 0; i < tdcdigihits.size(); i++) {
253 const DSCTDCDigiHit *digihit = tdcdigihits[i];
254
255 // Make sure sector is in valid range
256 if((digihit->sector <= 0) && (digihit->sector > MAX_SECTORS)) {
257 sprintf(str, "DSCDigiHit sector out of range! sector=%d (should be 1-%d)",
258 digihit->sector, MAX_SECTORS);
259 throw JException(str);
260 }
261
262 unsigned int id = digihit->sector - 1;
263 double T = locTTabUtilities->Convert_DigiTimeToNs_F1TDC(digihit) - tdc_time_offsets[id] + t_tdc_base;
264
265 // Look for existing hits to see if there is a match
266 // or create new one if there is no match
267 // Require that the trigger corrected TDC time fall within
268 // a reasonable time window so that when a hit is associated with
269 // a hit in the TDC and not the ADC it is a "decent" TDC hit
270 if (fabs(T) < HIT_TIME_WINDOW) {
271 DSCHit *hit = FindMatch(digihit->sector, T);
272 if (! hit) {
273 hit = new DSCHit;
274 hit->sector = digihit->sector;
275 hit->dE = 0.0;
276 hit->t_fADC= numeric_limits<double>::quiet_NaN();
277 hit->has_fADC=false;
278 _data.push_back(hit);
279 }
280
281 hit->has_TDC=true;
282 hit->t_TDC=T;
283
284 if (USE_TIMEWALK_CORRECTION && (hit->dE > 0.0) ) {
285 // Correct for time walk
286 // The correction is the form t=t_tdc- C1 (A^C2 - A0^C2)
287 // double A = hit->dE;
288 // double C1 = timewalk_parameters[id][1];
289 // double C2 = timewalk_parameters[id][2];
290 // double A0 = timewalk_parameters[id][3];
291 // T -= C1*(pow(A,C2) - pow(A0,C2));
292
293 // Correct for timewalk using pulse peak instead of pulse integral
294 double A = hit->pulse_height;
295 // double C0 = timewalk_parameters[id][0];
296 double C1 = timewalk_parameters[id][1];
297 double C2 = timewalk_parameters[id][2];
298 double A_THRESH = timewalk_parameters[id][3];
299 double A0 = timewalk_parameters[id][4];
300 // do correction
301 T -= C1*(pow(A/A_THRESH, C2) - pow(A0/A_THRESH, C2));
302 }
303
304 hit->t=T;
305
306 hit->AddAssociatedObject(digihit);
307 } // Hit time window cut
308
309 }
310
311
312 return NOERROR;
313}
314
315//------------------
316// FindMatch
317//------------------
318DSCHit* DSCHit_factory::FindMatch(int sector, double T)
319{
320 DSCHit *best_match=NULL__null;
321
322 // Loop over existing hits (from fADC) and look for a match
323 // in both the sector and the time.
324 for(unsigned int i = 0; i < _data.size(); i++)
325 {
326 DSCHit *hit = _data[i];
327
328 if (! isfinite(hit->t_fADC))
329 continue; // only match to fADC hits, not bachelor TDC hits
330
331 if (hit->sector != sector)
332 continue; // require identical sectors fired
333
334 double delta_T = fabs(hit->t - T);
335 if (delta_T > DELTA_T_ADC_TDC_MAX)
336 continue;
337
338 if (best_match != NULL__null)
339 {
340 if (delta_T < fabs(best_match->t - T))
341 best_match = hit;
342 } else best_match = hit;
343 }
344
345 return best_match;
346}
347
348//------------------
349// erun
350//------------------
351jerror_t DSCHit_factory::erun(void)
352{
353 return NOERROR;
354}
355
356//------------------
357// fini
358//------------------
359jerror_t DSCHit_factory::fini(void)
360{
361 return NOERROR;
362}
363
364
365//------------------------------------
366// GetConstant
367// Allow a few different interfaces
368//------------------------------------
369const double DSCHit_factory::GetConstant(const vector<double> &the_table,
370 const int in_sector) const
371{
372 char str[256];
373
374 if ( (in_sector < 0) || (in_sector >= MAX_SECTORS))
375 {
376 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
377 " requested=%d , should be %ud", in_sector, MAX_SECTORS);
378 cerr << str << endl;
379 throw JException(str);
380 }
381
382 return the_table[in_sector];
383}
384
385const double DSCHit_factory::GetConstant(const vector<double> &the_table,
386 const DSCDigiHit *in_digihit) const
387{
388 char str[256];
389
390 if ( (in_digihit->sector < 0) || (in_digihit->sector >= MAX_SECTORS))
391 {
392 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
393 " requested=%d , should be %ud",
394 in_digihit->sector, MAX_SECTORS);
395 cerr << str << endl;
396 throw JException(str);
397 }
398
399 return the_table[in_digihit->sector];
400}
401
402const double DSCHit_factory::GetConstant(const vector<double> &the_table,
403 const DSCHit *in_hit) const
404{
405
406 char str[256];
407
408 if ( (in_hit->sector < 0) || (in_hit->sector >= MAX_SECTORS))
409 {
410 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
411 " requested=%d , should be %ud",
412 in_hit->sector, MAX_SECTORS);
413 cerr << str << endl;
414 throw JException(str);
415 }
416
417 return the_table[in_hit->sector];
418}
419/*
420 const double DSCHit_factory::GetConstant(const vector<double> &the_table,
421 const DTranslationTable *ttab,
422 const int in_rocid,
423 const int in_slot,
424 const int in_channel) const
425 {
426 char str[256];
427
428 DTranslationTable::csc_t daq_index = { in_rocid, in_slot, in_channel };
429 DTranslationTable::DChannelInfo channel_info = ttab->GetDetectorIndex(daq_index);
430
431 if( (channel_info.sc.sector <= 0)
432 || (channel_info.sc.sector > static_cast<unsigned int>(MAX_SECTORS))) {
433 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
434 " requested=%d , should be %ud",
435 channel_info.sc.sector, MAX_SECTORS);
436 cerr << str << endl;
437 throw JException(str);
438 }
439
440 return the_table[channel_info.sc.sector];
441 }
442 */