Bug Summary

File:libraries/START_COUNTER/DSCHit_factory.cc
Location:line 202, 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 double ped_corr_pulse_peak = 0.;
173
174 // Make sure sector is in valid range
175 if( (digihit->sector <= 0) && (digihit->sector > MAX_SECTORS)) {
176 sprintf(str, "DSCDigiHit sector out of range! sector=%d (should be 1-%d)",
177 digihit->sector, MAX_SECTORS);
178 throw JException(str);
179 }
180
181 // Error checking for pre-Fall 2016 firmware
182 if(digihit->datasource == 1) {
183 // There is a slight difference between Mode 7 and 8 data
184 // The following condition signals an error state in the flash algorithm
185 // Do not make hits out of these
186 const Df250PulsePedestal* PPobj = NULL__null;
187 digihit->GetSingle(PPobj);
188 if (PPobj != NULL__null) {
189 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
190 } else
191 continue;
192
193 if (digihit->pulse_time == 0) continue; // Should already be caught, but I'll leave it
194 }
195
196 if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF))
197 continue;
198
199 // Initialize pedestal to one found in CCDB, but override it
200 // with one found in event if is available (?)
201 // For now, only keep events with a correct pedestal
202 double pedestal = a_pedestals[digihit->sector-1];
Value stored to 'pedestal' during its initialization is never read
203 double nsamples_integral = digihit->nsamples_integral;
204 double nsamples_pedestal = digihit->nsamples_pedestal;
205
206 // nsamples_pedestal should always be positive for valid data - err on the side of caution for now
207 if(nsamples_pedestal == 0) {
208 jerr << "DSCDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl;
209 continue;
210 }
211
212 if( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) {
213 pedestal = digihit->pedestal * (nsamples_integral/nsamples_pedestal);
214 } else {
215 continue;
216 }
217
218 double single_sample_ped = pedestal/nsamples_pedestal;
219
220 //if ( ((double)digihit->pulse_integral) < ADC_THRESHOLD) continue; // Will comment out until this is set to something useful by default
221
222 double A = (double)digihit->pulse_integral;
223 double T = (double)digihit->pulse_time;
224
225 DSCHit *hit = new DSCHit;
226 // Sectors are numbered from 1-30
227 hit->sector = digihit->sector;
228
229 hit->dE = a_scale * a_gains[digihit->sector-1] * (A - pedestal);
230 hit->t_fADC = t_scale * T - adc_time_offsets[hit->sector-1] + t_base;
231 hit->t_TDC = numeric_limits<double>::quiet_NaN();
232
233 hit->has_TDC = false;
234 hit->has_fADC = true;
235
236 hit->t = hit->t_fADC; // set time from fADC in case no TDC hit
237 hit->pulse_height = digihit->pulse_peak - single_sample_ped;
238
239 hit->AddAssociatedObject(digihit);
240
241 _data.push_back(hit);
242 }
243
244
245 // Next, loop over TDC hits, matching them to the
246 // existing fADC hits where possible and updating
247 // their time information. If no match is found, then
248 // create a new hit with just the TDC info.
249 vector<const DSCTDCDigiHit*> tdcdigihits;
250 loop->Get(tdcdigihits);
251 sort(tdcdigihits.begin(),tdcdigihits.end(),DSCHit_tdc_cmp);
252
253 for (unsigned int i = 0; i < tdcdigihits.size(); i++) {
254 const DSCTDCDigiHit *digihit = tdcdigihits[i];
255
256 // Make sure sector is in valid range
257 if((digihit->sector <= 0) && (digihit->sector > MAX_SECTORS)) {
258 sprintf(str, "DSCDigiHit sector out of range! sector=%d (should be 1-%d)",
259 digihit->sector, MAX_SECTORS);
260 throw JException(str);
261 }
262
263 unsigned int id = digihit->sector - 1;
264 double T = locTTabUtilities->Convert_DigiTimeToNs_F1TDC(digihit) - tdc_time_offsets[id] + t_tdc_base;
265
266 // Look for existing hits to see if there is a match
267 // or create new one if there is no match
268 // Require that the trigger corrected TDC time fall within
269 // a reasonable time window so that when a hit is associated with
270 // a hit in the TDC and not the ADC it is a "decent" TDC hit
271 if (fabs(T) < HIT_TIME_WINDOW) {
272 DSCHit *hit = FindMatch(digihit->sector, T);
273 if (! hit) {
274 hit = new DSCHit;
275 hit->sector = digihit->sector;
276 hit->dE = 0.0;
277 hit->t_fADC= numeric_limits<double>::quiet_NaN();
278 hit->has_fADC=false;
279 _data.push_back(hit);
280 }
281
282 hit->has_TDC=true;
283 hit->t_TDC=T;
284
285 if (USE_TIMEWALK_CORRECTION && (hit->dE > 0.0) ) {
286 // Correct for time walk
287 // The correction is the form t=t_tdc- C1 (A^C2 - A0^C2)
288 // double A = hit->dE;
289 // double C1 = timewalk_parameters[id][1];
290 // double C2 = timewalk_parameters[id][2];
291 // double A0 = timewalk_parameters[id][3];
292 // T -= C1*(pow(A,C2) - pow(A0,C2));
293
294 // Correct for timewalk using pulse peak instead of pulse integral
295 double A = hit->pulse_height;
296 // double C0 = timewalk_parameters[id][0];
297 double C1 = timewalk_parameters[id][1];
298 double C2 = timewalk_parameters[id][2];
299 double A_THRESH = timewalk_parameters[id][3];
300 double A0 = timewalk_parameters[id][4];
301 // do correction
302 T -= C1*(pow(A/A_THRESH, C2) - pow(A0/A_THRESH, C2));
303 }
304
305 hit->t=T;
306
307 hit->AddAssociatedObject(digihit);
308 } // Hit time window cut
309
310 }
311
312
313 return NOERROR;
314}
315
316//------------------
317// FindMatch
318//------------------
319DSCHit* DSCHit_factory::FindMatch(int sector, double T)
320{
321 DSCHit *best_match=NULL__null;
322
323 // Loop over existing hits (from fADC) and look for a match
324 // in both the sector and the time.
325 for(unsigned int i = 0; i < _data.size(); i++)
326 {
327 DSCHit *hit = _data[i];
328
329 if (! isfinite(hit->t_fADC))
330 continue; // only match to fADC hits, not bachelor TDC hits
331
332 if (hit->sector != sector)
333 continue; // require identical sectors fired
334
335 double delta_T = fabs(hit->t - T);
336 if (delta_T > DELTA_T_ADC_TDC_MAX)
337 continue;
338
339 if (best_match != NULL__null)
340 {
341 if (delta_T < fabs(best_match->t - T))
342 best_match = hit;
343 } else best_match = hit;
344 }
345
346 return best_match;
347}
348
349//------------------
350// erun
351//------------------
352jerror_t DSCHit_factory::erun(void)
353{
354 return NOERROR;
355}
356
357//------------------
358// fini
359//------------------
360jerror_t DSCHit_factory::fini(void)
361{
362 return NOERROR;
363}
364
365
366//------------------------------------
367// GetConstant
368// Allow a few different interfaces
369//------------------------------------
370const double DSCHit_factory::GetConstant(const vector<double> &the_table,
371 const int in_sector) const
372{
373 char str[256];
374
375 if ( (in_sector < 0) || (in_sector >= MAX_SECTORS))
376 {
377 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
378 " requested=%d , should be %ud", in_sector, MAX_SECTORS);
379 cerr << str << endl;
380 throw JException(str);
381 }
382
383 return the_table[in_sector];
384}
385
386const double DSCHit_factory::GetConstant(const vector<double> &the_table,
387 const DSCDigiHit *in_digihit) const
388{
389 char str[256];
390
391 if ( (in_digihit->sector < 0) || (in_digihit->sector >= MAX_SECTORS))
392 {
393 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
394 " requested=%d , should be %ud",
395 in_digihit->sector, MAX_SECTORS);
396 cerr << str << endl;
397 throw JException(str);
398 }
399
400 return the_table[in_digihit->sector];
401}
402
403const double DSCHit_factory::GetConstant(const vector<double> &the_table,
404 const DSCHit *in_hit) const
405{
406
407 char str[256];
408
409 if ( (in_hit->sector < 0) || (in_hit->sector >= MAX_SECTORS))
410 {
411 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
412 " requested=%d , should be %ud",
413 in_hit->sector, MAX_SECTORS);
414 cerr << str << endl;
415 throw JException(str);
416 }
417
418 return the_table[in_hit->sector];
419}
420/*
421 const double DSCHit_factory::GetConstant(const vector<double> &the_table,
422 const DTranslationTable *ttab,
423 const int in_rocid,
424 const int in_slot,
425 const int in_channel) const
426 {
427 char str[256];
428
429 DTranslationTable::csc_t daq_index = { in_rocid, in_slot, in_channel };
430 DTranslationTable::DChannelInfo channel_info = ttab->GetDetectorIndex(daq_index);
431
432 if( (channel_info.sc.sector <= 0)
433 || (channel_info.sc.sector > static_cast<unsigned int>(MAX_SECTORS))) {
434 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
435 " requested=%d , should be %ud",
436 channel_info.sc.sector, MAX_SECTORS);
437 cerr << str << endl;
438 throw JException(str);
439 }
440
441 return the_table[channel_info.sc.sector];
442 }
443 */