Bug Summary

File:libraries/TOF/DTOFHit_factory.cc
Location:line 203, column 16
Description:Value stored to 'pedestal' during its initialization is never read

Annotated Source Code

1// $Id$
2//
3// File: DTOFHit_factory.cc
4// Created: Wed Aug 7 09:30:17 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 <cmath>
12#include <vector>
13#include <limits>
14
15using namespace std;
16
17#include <TOF/DTOFDigiHit.h>
18#include <TOF/DTOFTDCDigiHit.h>
19#include "DTOFHit_factory.h"
20#include <DAQ/Df250PulseIntegral.h>
21#include <DAQ/Df250Config.h>
22#include <DAQ/DCODAROCInfo.h>
23using namespace jana;
24
25static bool COSMIC_DATA = false;
26
27//------------------
28// init
29//------------------
30jerror_t DTOFHit_factory::init(void)
31{
32 DELTA_T_ADC_TDC_MAX = 10.0; // ns
33 // DELTA_T_ADC_TDC_MAX = 30.0; // ns, value based on the studies from cosmic events
34 gPARMS->SetDefaultParameter("TOF:DELTA_T_ADC_TDC_MAX", DELTA_T_ADC_TDC_MAX,
35 "Maximum difference in ns between a (calibrated) fADC time and F1TDC time for them to be matched in a single hit");
36
37 int analyze_cosmic_data = 0;
38 gPARMS->SetDefaultParameter("TOF:COSMIC_DATA", analyze_cosmic_data,
39 "Special settings for analysing cosmic data");
40 if(analyze_cosmic_data > 0)
41 COSMIC_DATA = true;
42
43 CHECK_FADC_ERRORS = false;
44 gPARMS->SetDefaultParameter("TOF:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits");
45
46
47 /// Set basic conversion constants
48 a_scale = 0.2/5.2E5;
49 t_scale = 0.0625; // 62.5 ps/count
50 t_base = 0.; // ns
51 t_base_tdc = 0.; // ns
52
53 if(COSMIC_DATA)
54 // Hardcoding of 110 taken from cosmics events
55 tdc_adc_time_offset = 110.;
56 else
57 tdc_adc_time_offset = 0.;
58
59 TOF_NUM_PLANES = 2;
60 TOF_NUM_BARS = 44;
61
62 return NOERROR;
63}
64
65//------------------
66// brun
67//------------------
68jerror_t DTOFHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
69{
70 // Only print messages for one thread whenever run number change
71 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
72 static set<int> runs_announced;
73 pthread_mutex_lock(&print_mutex);
74 bool print_messages = false;
75 if(runs_announced.find(runnumber) == runs_announced.end()){
76 print_messages = true;
77 runs_announced.insert(runnumber);
78 }
79 pthread_mutex_unlock(&print_mutex);
80
81 // read in geometry information
82 vector<const DTOFGeometry*> tofGeomVect;
83 eventLoop->Get( tofGeomVect );
84 if(tofGeomVect.size()<1) return OBJECT_NOT_AVAILABLE;
85 const DTOFGeometry& tofGeom = *(tofGeomVect[0]);
86
87 /// Read in calibration constants
88 vector<double> raw_adc_pedestals;
89 vector<double> raw_adc_gains;
90 vector<double> raw_adc_offsets;
91 vector<double> raw_tdc_offsets;
92 vector<double> raw_adc2E;
93
94 if(print_messages) jout << "In DTOFHit_factory, loading constants..." << endl;
95
96 // load scale factors
97 map<string,double> scale_factors;
98 if(eventLoop->GetCalib("/TOF/digi_scales", scale_factors))
99 jout << "Error loading /TOF/digi_scales !" << endl;
100 if( scale_factors.find("TOF_ADC_ASCALE") != scale_factors.end() ) {
101 ; //a_scale = scale_factors["TOF_ADC_ASCALE"];
102 } else {
103 jerr << "Unable to get TOF_ADC_ASCALE from /TOF/digi_scales !" << endl;
104 }
105 if( scale_factors.find("TOF_ADC_TSCALE") != scale_factors.end() ) {
106 ; //t_scale = scale_factors["TOF_ADC_TSCALE"];
107 } else {
108 jerr << "Unable to get TOF_ADC_TSCALE from /TOF/digi_scales !" << endl;
109 }
110
111 // load base time offset
112 map<string,double> base_time_offset;
113 if (eventLoop->GetCalib("/TOF/base_time_offset",base_time_offset))
114 jout << "Error loading /TOF/base_time_offset !" << endl;
115 if (base_time_offset.find("TOF_BASE_TIME_OFFSET") != base_time_offset.end())
116 t_base = base_time_offset["TOF_BASE_TIME_OFFSET"];
117 else
118 jerr << "Unable to get TOF_BASE_TIME_OFFSET from /TOF/base_time_offset !" << endl;
119
120 if (base_time_offset.find("TOF_TDC_BASE_TIME_OFFSET") != base_time_offset.end())
121 t_base_tdc = base_time_offset["TOF_TDC_BASE_TIME_OFFSET"];
122 else
123 jerr << "Unable to get TOF_TDC_BASE_TIME_OFFSET from /TOF/base_time_offset !" << endl;
124
125 // load constant tables
126 if(eventLoop->GetCalib("TOF/pedestals", raw_adc_pedestals))
127 jout << "Error loading /TOF/pedestals !" << endl;
128 if(eventLoop->GetCalib("TOF/gains", raw_adc_gains))
129 jout << "Error loading /TOF/gains !" << endl;
130 if(eventLoop->GetCalib("TOF/adc_timing_offsets", raw_adc_offsets))
131 jout << "Error loading /TOF/adc_timing_offsets !" << endl;
132 if(eventLoop->GetCalib("TOF/timing_offsets", raw_tdc_offsets))
133 jout << "Error loading /TOF/timing_offsets !" << endl;
134 if(eventLoop->GetCalib("TOF/timewalk_parms", timewalk_parameters))
135 jout << "Error loading /TOF/timewalk_parms !" << endl;
136
137 FillCalibTable(adc_pedestals, raw_adc_pedestals, tofGeom);
138 FillCalibTable(adc_gains, raw_adc_gains, tofGeom);
139 FillCalibTable(adc_time_offsets, raw_adc_offsets, tofGeom);
140 FillCalibTable(tdc_time_offsets, raw_tdc_offsets, tofGeom);
141
142 if(eventLoop->GetCalib("TOF/adc2E", raw_adc2E))
143 jout << "Error loading /TOF/adc2E !" << endl;
144
145 for (unsigned int n=0; n<raw_adc2E.size(); n++){
146 adc2E[n] = raw_adc2E[n];
147 }
148
149 /*
150 CheckCalibTable(adc_pedestals,"/TOF/pedestals");
151 CheckCalibTable(adc_gains,"/TOF/gains");
152 CheckCalibTable(adc_time_offsets,"/TOF/adc_timing_offsets");
153 CheckCalibTable(tdc_time_offsets,"/TOF/timing_offsets");
154 */
155
156 return NOERROR;
157}
158
159//------------------
160// evnt
161//------------------
162jerror_t DTOFHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
163{
164 /// Generate DTOFHit object for each DTOFDigiHit object.
165 /// This is where the first set of calibration constants
166 /// is applied to convert from digitzed units into natural
167 /// units.
168 ///
169 /// Note that this code does NOT get called for simulated
170 /// data in HDDM format. The HDDM event source will copy
171 /// the precalibrated values directly into the _data vector.
172
173 const DTTabUtilities* locTTabUtilities = NULL__null;
174 loop->GetSingle(locTTabUtilities);
175
176 // First, make hits out of all fADC250 hits
177 vector<const DTOFDigiHit*> digihits;
178 loop->Get(digihits);
179 for(unsigned int i=0; i<digihits.size(); i++){
180 const DTOFDigiHit *digihit = digihits[i];
181
182 // Error checking for pre-Fall 2016 firmware
183 if(digihit->datasource == 1) {
184 // There is a slight difference between Mode 7 and 8 data
185 // The following condition signals an error state in the flash algorithm
186 // Do not make hits out of these
187 const Df250PulsePedestal* PPobj = NULL__null;
188 digihit->GetSingle(PPobj);
189 if (PPobj != NULL__null) {
190 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
191 } else
192 continue;
193
194 //if (digihit->pulse_time == 0) continue; // Should already be caught
195 }
196
197 if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF))
198 continue;
199
200 // Initialize pedestal to one found in CCDB, but override it
201 // with one found in event if is available (?)
202 // For now, only keep events with a correct pedestal
203 double pedestal = GetConstant(adc_pedestals, digihit);
Value stored to 'pedestal' during its initialization is never read
204 double nsamples_integral = digihit->nsamples_integral;
205 double nsamples_pedestal = digihit->nsamples_pedestal;
206
207 // nsamples_pedestal should always be positive for valid data - err on the side of caution for now
208 if(nsamples_pedestal == 0) {
209 jerr << "DTOFDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl;
210 continue;
211 }
212
213 if( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) {
214 pedestal = digihit->pedestal * (nsamples_integral/nsamples_pedestal);
215 } else {
216 continue;
217 }
218
219 //double single_sample_ped = pedestal/nsamples_pedestal;
220
221 // Apply calibration constants here
222 double A = (double)digihit->pulse_integral;
223 double T = (double)digihit->pulse_time;
224 T = t_scale * T - GetConstant(adc_time_offsets, digihit) + t_base;
225 double dA = A - pedestal;
226
227 if (dA<0) continue;
228
229 DTOFHit *hit = new DTOFHit;
230 hit->plane = digihit->plane;
231 hit->bar = digihit->bar;
232 hit->end = digihit->end;
233 hit->dE=dA; // this will be scaled to energy units later
234
235 if(COSMIC_DATA)
236 hit->dE = (A - 55*pedestal); // value of 55 is taken from (NSB,NSA)=(10,45) in the confg file
237
238 hit->t_TDC=numeric_limits<double>::quiet_NaN();
239 hit->t_fADC=T;
240 hit->t = hit->t_fADC; // set initial time to the ADC time, in case there's no matching TDC hit
241
242 hit->has_fADC=true;
243 hit->has_TDC=false;
244
245 /*
246 cout << "TOF ADC hit = (" << hit->plane << "," << hit->bar << "," << hit->end << ") "
247 << t_scale << " " << T << " "
248 << GetConstant(adc_time_offsets, digihit) << " "
249 << t_scale*GetConstant(adc_time_offsets, digihit) << " " << hit->t << endl;
250 */
251
252 hit->AddAssociatedObject(digihit);
253
254 _data.push_back(hit);
255 }
256
257 //Get the TDC hits
258 vector<const DTOFTDCDigiHit*> tdcdigihits;
259 loop->Get(tdcdigihits);
260
261 // Next, loop over TDC hits, matching them to the
262 // existing fADC hits where possible and updating
263 // their time information. If no match is found, then
264 // create a new hit with just the TDC info.
265 for(unsigned int i=0; i<tdcdigihits.size(); i++)
266 {
267 const DTOFTDCDigiHit *digihit = tdcdigihits[i];
268
269 // Apply calibration constants here
270 double T = locTTabUtilities->Convert_DigiTimeToNs_CAEN1290TDC(digihit);
271 T += t_base_tdc - GetConstant(tdc_time_offsets, digihit) + tdc_adc_time_offset;
272
273 /*
274 cout << "TOF TDC hit = (" << digihit->plane << "," << digihit->bar << "," << digihit->end << ") "
275 << T << " " << GetConstant(tdc_time_offsets, digihit) << endl;
276 */
277
278 // Look for existing hits to see if there is a match
279 // or create new one if there is no match
280 DTOFHit *hit = FindMatch(digihit->plane, digihit->bar, digihit->end, T);
281 //DTOFHit *hit = FindMatch(digihit->plane, hit->bar, hit->end, T);
282 if(!hit){
283 hit = new DTOFHit;
284 hit->plane = digihit->plane;
285 hit->bar = digihit->bar;
286 hit->end = digihit->end;
287 hit->dE = 0.0;
288 hit->t_fADC=numeric_limits<double>::quiet_NaN();
289 hit->has_fADC=false;
290
291 _data.push_back(hit);
292 }
293 hit->has_TDC=true;
294 hit->t_TDC=T;
295
296 if (hit->dE>0.){
297 // time walk correction
298 // The correction is the form t=t_tdc- C1 (A^C2 - A0^C2)
299 int id=88*hit->plane+44*hit->end+hit->bar-1;
300 double A=hit->dE;
301 double C1=timewalk_parameters[id][1];
302 double C2=timewalk_parameters[id][2];
303 double A0=timewalk_parameters[id][3];
304 T-=C1*(pow(A,C2)-pow(A0,C2));
305 }
306 hit->t=T;
307
308 hit->AddAssociatedObject(digihit);
309 }
310
311 // Apply calibration constants to convert pulse integrals to energy units
312 for (unsigned int i=0;i<_data.size();i++){
313 int id=88*_data[i]->plane + 44*_data[i]->end + _data[i]->bar-1;
314 _data[i]->dE *= adc2E[id];
315 //cout<<id<<" "<< adc2E[id]<<" "<<_data[i]->dE<<endl;
316 }
317
318 return NOERROR;
319}
320
321//------------------
322// FindMatch
323//------------------
324DTOFHit* DTOFHit_factory::FindMatch(int plane, int bar, int end, double T)
325{
326 DTOFHit* best_match = NULL__null;
327
328 // Loop over existing hits (from fADC) and look for a match
329 // in both the sector and the time.
330 for(unsigned int i=0; i<_data.size(); i++){
331 DTOFHit *hit = _data[i];
332
333 if(!isfinite(hit->t_fADC)) continue; // only match to fADC hits, not bachelor TDC hits
334 if(hit->plane != plane) continue;
335 if(hit->bar != bar) continue;
336 if(hit->end != end) continue;
337
338 //double delta_T = fabs(hit->t - T);
339 double delta_T = fabs(T - hit->t);
340 if(delta_T > DELTA_T_ADC_TDC_MAX) continue;
341
342 return hit;
343
344 // if there are multiple hits, pick the one that is closest in time
345 if(best_match != NULL__null) {
346 if(delta_T < fabs(best_match->t - T))
347 best_match = hit;
348 } else {
349 best_match = hit;
350 }
351
352 }
353
354 return best_match;
355}
356
357//------------------
358// erun
359//------------------
360jerror_t DTOFHit_factory::erun(void)
361{
362 return NOERROR;
363}
364
365//------------------
366// fini
367//------------------
368jerror_t DTOFHit_factory::fini(void)
369{
370 return NOERROR;
371}
372
373
374//------------------
375// FillCalibTable
376//------------------
377void DTOFHit_factory::FillCalibTable(tof_digi_constants_t &table, vector<double> &raw_table,
378 const DTOFGeometry &tofGeom)
379{
380 char str[256];
381 int channel = 0;
382
383 // reset the table before filling it
384 table.clear();
385
386 for(int plane=0; plane<tofGeom.NLAYERS; plane++) {
387 int plane_index=2*TOF_NUM_BARS*plane;
388 table.push_back( vector< pair<double,double> >(TOF_NUM_BARS) );
389 for(int bar=0; bar<TOF_NUM_BARS; bar++) {
390 table[plane][bar]
391 = pair<double,double>(raw_table[plane_index+bar],
392 raw_table[plane_index+TOF_NUM_BARS+bar]);
393 channel+=2;
394 }
395 }
396
397 // check to make sure that we loaded enough channels
398 if(channel != TOF_MAX_CHANNELS) {
399 sprintf(str, "Not enough channels for TOF table! channel=%d (should be %d)",
400 channel, TOF_MAX_CHANNELS);
401 cerr << str << endl;
402 throw JException(str);
403 }
404}
405
406
407//------------------------------------
408// GetConstant
409// Allow a few different interfaces
410// NOTE: LoadGeometry() must be called before calling these functions
411//
412// TOF Geometry as defined in the Translation Table:
413// plane = 0-1
414// bar = 1-44
415// end = 0-1
416// Note the different counting schemes used
417//------------------------------------
418const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table,
419 const int in_plane, const int in_bar, const int in_end) const
420{
421 char str[256];
422
423 if( (in_plane < 0) || (in_plane > TOF_NUM_PLANES)) {
424 sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_plane, TOF_NUM_PLANES);
425 cerr << str << endl;
426 throw JException(str);
427 }
428 if( (in_bar <= 0) || (in_bar > TOF_NUM_BARS)) {
429 sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_bar, TOF_NUM_BARS);
430 cerr << str << endl;
431 throw JException(str);
432 }
433 if( (in_end != 0) && (in_end != 1) ) {
434 sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_end);
435 cerr << str << endl;
436 throw JException(str);
437 }
438
439 // we have two ends, indexed as 0/1
440 // could be north/south or up/down depending on the bar orientation
441 if(in_end == 0) {
442 return the_table[in_plane][in_bar].first;
443 } else {
444 return the_table[in_plane][in_bar].second;
445 }
446}
447
448const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table,
449 const DTOFHit *in_hit) const
450{
451 char str[256];
452
453 if( (in_hit->plane < 0) || (in_hit->plane > TOF_NUM_PLANES)) {
454 sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_hit->plane, TOF_NUM_PLANES);
455 cerr << str << endl;
456 throw JException(str);
457 }
458 if( (in_hit->bar <= 0) || (in_hit->bar > TOF_NUM_BARS)) {
459 sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_hit->bar, TOF_NUM_BARS);
460 cerr << str << endl;
461 throw JException(str);
462 }
463 if( (in_hit->end != 0) && (in_hit->end != 1) ) {
464 sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_hit->end);
465 cerr << str << endl;
466 throw JException(str);
467 }
468
469 // we have two ends, indexed as 0/1
470 // could be north/south or up/down depending on the bar orientation
471 if(in_hit->end == 0) {
472 return the_table[in_hit->plane][in_hit->bar-1].first;
473 } else {
474 return the_table[in_hit->plane][in_hit->bar-1].second;
475 }
476}
477
478const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table,
479 const DTOFDigiHit *in_digihit) const
480{
481 char str[256];
482
483 if( (in_digihit->plane < 0) || (in_digihit->plane > TOF_NUM_PLANES)) {
484 sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->plane, TOF_NUM_PLANES);
485 cerr << str << endl;
486 throw JException(str);
487 }
488 if( (in_digihit->bar <= 0) || (in_digihit->bar > TOF_NUM_BARS)) {
489 sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->bar, TOF_NUM_BARS);
490 cerr << str << endl;
491 throw JException(str);
492 }
493 if( (in_digihit->end != 0) && (in_digihit->end != 1) ) {
494 sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_digihit->end);
495 cerr << str << endl;
496 throw JException(str);
497 }
498
499 // we have two ends, indexed as 0/1
500 // could be north/south or up/down depending on the bar orientation
501 if(in_digihit->end == 0) {
502 return the_table[in_digihit->plane][in_digihit->bar-1].first;
503 } else {
504 return the_table[in_digihit->plane][in_digihit->bar-1].second;
505 }
506}
507
508const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table,
509 const DTOFTDCDigiHit *in_digihit) const
510{
511 char str[256];
512
513 if( (in_digihit->plane < 0) || (in_digihit->plane > TOF_NUM_PLANES)) {
514 sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->plane, TOF_NUM_PLANES);
515 cerr << str << endl;
516 throw JException(str);
517 }
518 if( (in_digihit->bar <= 0) || (in_digihit->bar > TOF_NUM_BARS)) {
519 sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->bar, TOF_NUM_BARS);
520 cerr << str << endl;
521 throw JException(str);
522 }
523 if( (in_digihit->end != 0) && (in_digihit->end != 1) ) {
524 sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_digihit->end);
525 cerr << str << endl;
526 throw JException(str);
527 }
528
529 // we have two ends, indexed as 0/1
530 // could be north/south or up/down depending on the bar orientation
531 if(in_digihit->end == 0) {
532 return the_table[in_digihit->plane][in_digihit->bar-1].first;
533 } else {
534 return the_table[in_digihit->plane][in_digihit->bar-1].second;
535 }
536}
537
538/*
539 const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table,
540 const DTranslationTable *ttab,
541 const int in_rocid, const int in_slot, const int in_channel) const {
542
543 char str[256];
544
545 DTranslationTable::csc_t daq_index = { in_rocid, in_slot, in_channel };
546 DTranslationTable::DChannelInfo channel_info = ttab->GetDetectorIndex(daq_index);
547
548 if( (channel_info.tof.plane <= 0)
549 || (channel_info.tof.plane > static_cast<unsigned int>(TOF_NUM_PLANES))) {
550 sprintf(str, "Bad plane # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", channel_info.tof.plane, TOF_NUM_PLANES);
551 cerr << str << endl;
552 throw JException(str);
553 }
554 if( (channel_info.tof.bar <= 0)
555 || (channel_info.tof.bar > static_cast<unsigned int>(TOF_NUM_BARS))) {
556 sprintf(str, "Bad bar # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", channel_info.tof.bar, TOF_NUM_BARS);
557 cerr << str << endl;
558 throw JException(str);
559 }
560 if( (channel_info.tof.end != 0) && (channel_info.tof.end != 1) ) {
561 sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", channel_info.tof.end);
562 cerr << str << endl;
563 throw JException(str);
564 }
565
566 int the_cell = DTOFGeometry::cellId(channel_info.tof.module,
567 channel_info.tof.layer,
568 channel_info.tof.sector);
569
570 if(channel_info.tof.end == DTOFGeometry::kUpstream) {
571// handle the upstream end
572return the_table.at(the_cell).first;
573} else {
574// handle the downstream end
575return the_table.at(the_cell).second;
576}
577}
578*/