Bug Summary

File:libraries/START_COUNTER/DSCHit_factory.cc
Location:line 215, 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 return NOERROR;
69}
70
71//------------------
72// brun
73//------------------
74jerror_t DSCHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
75{
76 // Only print messages for one thread whenever run number change
77 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
78 static set<int> runs_announced;
79 pthread_mutex_lock(&print_mutex);
80 bool print_messages = false;
81 if(runs_announced.find(runnumber) == runs_announced.end()){
82 print_messages = true;
83 runs_announced.insert(runnumber);
84 }
85 pthread_mutex_unlock(&print_mutex);
86
87 /// Read in calibration constants
88 if(print_messages) jout << "In DSCHit_factory, loading constants..." << endl;
89
90 // load scale factors
91 map<string,double> scale_factors;
92 // a_scale (SC_ADC_SCALE)
93 if (eventLoop->GetCalib("/START_COUNTER/digi_scales", scale_factors))
94 jout << "Error loading /START_COUNTER/digi_scales !" << endl;
95 if (scale_factors.find("SC_ADC_ASCALE") != scale_factors.end())
96 a_scale = scale_factors["SC_ADC_ASCALE"];
97 else
98 jerr << "Unable to get SC_ADC_ASCALE from /START_COUNTER/digi_scales !"
99 << endl;
100 // t_scale (SC_ADC_SCALE)
101 if (scale_factors.find("SC_ADC_TSCALE") != scale_factors.end())
102 t_scale = scale_factors["SC_ADC_TSCALE"];
103 else
104 jerr << "Unable to get SC_ADC_TSCALE from /START_COUNTER/digi_scales !"
105 << endl;
106
107 // load base time offset
108 map<string,double> base_time_offset;
109 // t_base (SC_BASE_TIME_OFFSET)
110 if (eventLoop->GetCalib("/START_COUNTER/base_time_offset",base_time_offset))
111 jout << "Error loading /START_COUNTER/base_time_offset !" << endl;
112 if (base_time_offset.find("SC_BASE_TIME_OFFSET") != base_time_offset.end())
113 t_base = base_time_offset["SC_BASE_TIME_OFFSET"];
114 else
115 jerr << "Unable to get SC_BASE_TIME_OFFSET from /START_COUNTER/base_time_offset !" << endl;
116 // t_tdc_base (SC_TDC_BASE_TIME_OFFSET)
117 if (base_time_offset.find("SC_TDC_BASE_TIME_OFFSET") != base_time_offset.end())
118 t_tdc_base = base_time_offset["SC_TDC_BASE_TIME_OFFSET"];
119 else
120 jerr << "Unable to get SC_BASE_TIME_OFFSET from /START_COUNTER/base_time_offset !" << endl;
121
122 // load constant tables
123 // a_gains (gains)
124 if (eventLoop->GetCalib("/START_COUNTER/gains", a_gains))
125 jout << "Error loading /START_COUNTER/gains !" << endl;
126 // a_pedestals (pedestals)
127 if (eventLoop->GetCalib("/START_COUNTER/pedestals", a_pedestals))
128 jout << "Error loading /START_COUNTER/pedestals !" << endl;
129 // adc_time_offsets (adc_timing_offsets)
130 if (eventLoop->GetCalib("/START_COUNTER/adc_timing_offsets", adc_time_offsets))
131 jout << "Error loading /START_COUNTER/adc_timing_offsets !" << endl;
132 // tdc_time_offsets (tdc_timing_offsets)
133 if (eventLoop->GetCalib("/START_COUNTER/tdc_timing_offsets", tdc_time_offsets))
134 jout << "Error loading /START_COUNTER/tdc_timing_offsets !" << endl;
135 // timewalk_parameters (timewalk_parms)
136 if(eventLoop->GetCalib("START_COUNTER/timewalk_parms_v2", timewalk_parameters))
137 jout << "Error loading /START_COUNTER/timewalk_parms_v2 !" << endl;
138
139
140 /*
141 // load higher order corrections
142 map<string,double> in_prop_corr;
143 map<string,double> in_atten_corr;
144
145 if(!eventLoop->GetCalib("/START_COUNTER/propagation_speed",in_prop_corr ))
146 jout << "Error loading /START_COUNTER/propagation_speed !" << endl;
147 if(!eventLoop->GetCalib("/START_COUNTER/attenuation_factor", in_atten_corr))
148 jout << "Error loading /START_COUNTER/attenuation_factor !" << endl;
149
150 // propogation correction: A + Bx
151 propogation_corr_factors.push_back(in_prop_corr["A"]);
152 propogation_corr_factors.push_back(in_prop_corr["B"]);
153
154 // attenuation correction: A + Bx + Cx^2 + Dx^3 + Ex^4 + Fx^5
155 attenuation_corr_factors.push_back(in_atten_corr["A"]);
156 attenuation_corr_factors.push_back(in_atten_corr["B"]);
157 attenuation_corr_factors.push_back(in_atten_corr["C"]);
158 attenuation_corr_factors.push_back(in_atten_corr["D"]);
159 attenuation_corr_factors.push_back(in_atten_corr["E"]);
160 attenuation_corr_factors.push_back(in_atten_corr["F"]);
161 */
162
163 return NOERROR;
164}
165
166//------------------
167// evnt
168//------------------
169jerror_t DSCHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
170{
171 /// Generate DSCHit object for each DSCDigiHit object.
172 /// This is where the first set of calibration constants
173 /// is applied to convert from digitzed units into natural
174 /// units.
175 ///
176 /// Note that this code does NOT get called for simulated
177 /// data in HDDM format. The HDDM event source will copy
178 /// the precalibrated values directly into the _data vector.
179
180 // First, make hits out of all fADC250 hits
181 vector<const DSCDigiHit*> digihits;
182 loop->Get(digihits);
183 sort(digihits.begin(),digihits.end(),DSCHit_fadc_cmp);
184
185 const DTTabUtilities* locTTabUtilities = NULL__null;
186 loop->GetSingle(locTTabUtilities);
187
188 char str[256];
189
190 for (unsigned int i = 0; i < digihits.size(); i++)
191 {
192 const DSCDigiHit *digihit = digihits[i];
193
194 // There is a slight difference between Mode 7 and 8 data
195 // The following condition signals an error state in the flash algorithm
196 // Do not make hits out of these
197 const Df250PulsePedestal* PPobj = NULL__null;
198 digihit->GetSingle(PPobj);
199 if (PPobj != NULL__null)
200 {
201 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
202 }
203 else continue;
204
205 // Make sure sector is in valid range
206 if( (digihit->sector <= 0) && (digihit->sector > MAX_SECTORS))
207 {
208 sprintf(str, "DSCDigiHit sector out of range! sector=%d (should be 1-%d)",
209 digihit->sector, MAX_SECTORS);
210 throw JException(str);
211 }
212
213 // Initialize pedestal to one found in CCDB, but override it
214 // with one found in event if is available
215 double pedestal = a_pedestals[digihit->sector-1];
Value stored to 'pedestal' during its initialization is never read
216 const Df250PulseIntegral *pulse_integral = NULL__null;
217 digihit->GetSingle(pulse_integral);
218 if (pulse_integral != NULL__null)
219 {
220 double single_sample_ped = (double)pulse_integral->pedestal;
221 double nsamples_integral = (double)pulse_integral->nsamples_integral;
222 double nsamples_pedestal = (double)pulse_integral->nsamples_pedestal;
223 pedestal = single_sample_ped * nsamples_integral/nsamples_pedestal;
224 }
225 else continue;
226
227 // Apply calibration constants here
228 double A = (double)digihit->pulse_integral;
229 double T = (double)digihit->pulse_time;
230 double dA = A - pedestal;
231 double ped_corr_pulse_peak = PPobj->pulse_peak - PPobj->pedestal;
232
233 if (digihit->pulse_time == 0) continue; // Should already be caught, but I'll leave it
234 //if (dA < ADC_THRESHOLD) continue; // Will comment out until this is set to something useful by default
235
236 DSCHit *hit = new DSCHit;
237 hit->sector = digihit->sector;
238
239 // Sectors are numbered from 1-30
240 hit->dE = dA; // This will be scaled to energy units later
241 hit->t_fADC = t_scale * T - adc_time_offsets[hit->sector-1] + t_base;
242 hit->t_TDC = numeric_limits<double>::quiet_NaN();
243
244 hit->has_TDC = false;
245 hit->has_fADC = true;
246
247 hit->t = hit->t_fADC; // set time from fADC in case no TDC hit
248 hit->pulse_height = ped_corr_pulse_peak;
249
250 // add in higher order corrections?
251
252 hit->AddAssociatedObject(digihit);
253
254 _data.push_back(hit);
255 }
256
257
258 // Get the trigger time from the f1 TDC
259 vector<const DF1TDCHit*> tdchit;
260 eventLoop->Get(tdchit);
261
262 // Next, loop over TDC hits, matching them to the
263 // existing fADC hits where possible and updating
264 // their time information. If no match is found, then
265 // create a new hit with just the TDC info.
266 vector<const DSCTDCDigiHit*> tdcdigihits;
267 loop->Get(tdcdigihits);
268 sort(tdcdigihits.begin(),tdcdigihits.end(),DSCHit_tdc_cmp);
269
270 for (unsigned int i = 0; i < tdcdigihits.size(); i++)
271 {
272 const DSCTDCDigiHit *digihit = tdcdigihits[i];
273
274 // Make sure sector is in valid range
275 if((digihit->sector <= 0) && (digihit->sector > MAX_SECTORS))
276 {
277 sprintf(str, "DSCDigiHit sector out of range! sector=%d (should be 1-%d)",
278 digihit->sector, MAX_SECTORS);
279 throw JException(str);
280 }
281
282 unsigned int id = digihit->sector - 1;
283 double T = locTTabUtilities->Convert_DigiTimeToNs_F1TDC(digihit) - tdc_time_offsets[id] + t_tdc_base;
284
285 // cout << "T = " << T << endl;
286 // jout << "T = " << T << endl;
287
288 // Look for existing hits to see if there is a match
289 // or create new one if there is no match
290 // Require that the trigger corrected TDC time fall within
291 // a reasonable time window so that when a hit is associated with
292 // a hit in the TDC and not the ADC it is a "decent" TDC hit
293 if (fabs(T) < HIT_TIME_WINDOW)
294 {
295 //jout << " T cut = " << T << endl;
296 DSCHit *hit = FindMatch(digihit->sector, T);
297 if (! hit)
298 {
299 hit = new DSCHit;
300 hit->sector = digihit->sector;
301 hit->dE = 0.0;
302 hit->t_fADC= numeric_limits<double>::quiet_NaN();
303 hit->has_fADC=false;
304 _data.push_back(hit);
305 }
306
307 hit->has_TDC=true;
308 hit->t_TDC=T;
309 //jout << "t_tDC = " << hit->t_TDC << endl;
310
311
312 if (USE_TIMEWALK_CORRECTION && (hit->dE > 0.0) )
313 {
314 // // Correct for time walk
315 // // The correction is the form t=t_tdc- C1 (A^C2 - A0^C2)
316 // double A = hit->dE;
317 // double C1 = timewalk_parameters[id][1];
318 // double C2 = timewalk_parameters[id][2];
319 // double A0 = timewalk_parameters[id][3];
320 // T -= C1*(pow(A,C2) - pow(A0,C2));
321
322 // Correct for timewalk using pulse peak instead of pulse integral
323 double A = hit->pulse_height;
324 // double C0 = timewalk_parameters[id][0];
325 double C1 = timewalk_parameters[id][1];
326 double C2 = timewalk_parameters[id][2];
327 double A_THRESH = timewalk_parameters[id][3];
328 double A0 = timewalk_parameters[id][4];
329 // do correction
330 T -= C1*(pow(A/A_THRESH, C2) - pow(A0/A_THRESH, C2));
331 }
332 hit->t=T;
333 //jout << " T cut TW Corr = " << T << endl;
334 hit->AddAssociatedObject(digihit);
335 } // Hit time window cut
336
337 }
338
339
340 // Apply calibration constants to convert pulse integrals to energy
341 // units
342 for (unsigned int i=0;i<_data.size();i++)
343 {
344 _data[i]->dE*=a_scale * a_gains[_data[i]->sector-1];
345 }
346
347
348 return NOERROR;
349}
350
351//------------------
352// FindMatch
353//------------------
354DSCHit* DSCHit_factory::FindMatch(int sector, double T)
355{
356 DSCHit *best_match=NULL__null;
357
358 // Loop over existing hits (from fADC) and look for a match
359 // in both the sector and the time.
360 for(unsigned int i = 0; i < _data.size(); i++)
361 {
362 DSCHit *hit = _data[i];
363
364 if (! isfinite(hit->t_fADC))
365 continue; // only match to fADC hits, not bachelor TDC hits
366
367 if (hit->sector != sector)
368 continue; // require identical sectors fired
369
370 double delta_T = fabs(hit->t - T);
371 if (delta_T > DELTA_T_ADC_TDC_MAX)
372 continue;
373
374 if (best_match != NULL__null)
375 {
376 if (delta_T < fabs(best_match->t - T))
377 best_match = hit;
378 } else best_match = hit;
379 }
380
381 return best_match;
382}
383
384//------------------
385// erun
386//------------------
387jerror_t DSCHit_factory::erun(void)
388{
389 return NOERROR;
390}
391
392//------------------
393// fini
394//------------------
395jerror_t DSCHit_factory::fini(void)
396{
397 return NOERROR;
398}
399
400
401//------------------------------------
402// GetConstant
403// Allow a few different interfaces
404//------------------------------------
405const double DSCHit_factory::GetConstant(const vector<double> &the_table,
406 const int in_sector) const
407{
408 char str[256];
409
410 if ( (in_sector < 0) || (in_sector >= MAX_SECTORS))
411 {
412 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
413 " requested=%d , should be %ud", in_sector, MAX_SECTORS);
414 cerr << str << endl;
415 throw JException(str);
416 }
417
418 return the_table[in_sector];
419}
420
421const double DSCHit_factory::GetConstant(const vector<double> &the_table,
422 const DSCDigiHit *in_digihit) const
423{
424 char str[256];
425
426 if ( (in_digihit->sector < 0) || (in_digihit->sector >= MAX_SECTORS))
427 {
428 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
429 " requested=%d , should be %ud",
430 in_digihit->sector, MAX_SECTORS);
431 cerr << str << endl;
432 throw JException(str);
433 }
434
435 return the_table[in_digihit->sector];
436}
437
438const double DSCHit_factory::GetConstant(const vector<double> &the_table,
439 const DSCHit *in_hit) const
440{
441
442 char str[256];
443
444 if ( (in_hit->sector < 0) || (in_hit->sector >= MAX_SECTORS))
445 {
446 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
447 " requested=%d , should be %ud",
448 in_hit->sector, MAX_SECTORS);
449 cerr << str << endl;
450 throw JException(str);
451 }
452
453 return the_table[in_hit->sector];
454}
455/*
456 const double DSCHit_factory::GetConstant(const vector<double> &the_table,
457 const DTranslationTable *ttab,
458 const int in_rocid,
459 const int in_slot,
460 const int in_channel) const
461 {
462 char str[256];
463
464 DTranslationTable::csc_t daq_index = { in_rocid, in_slot, in_channel };
465 DTranslationTable::DChannelInfo channel_info = ttab->GetDetectorIndex(daq_index);
466
467 if( (channel_info.sc.sector <= 0)
468 || (channel_info.sc.sector > static_cast<unsigned int>(MAX_SECTORS))) {
469 sprintf(str, "Bad sector # requested in DSCHit_factory::GetConstant()!"
470 " requested=%d , should be %ud",
471 channel_info.sc.sector, MAX_SECTORS);
472 cerr << str << endl;
473 throw JException(str);
474 }
475
476 return the_table[channel_info.sc.sector];
477 }
478 */