Bug Summary

File:libraries/CDC/DCDCHit_factory.cc
Location:line 207, column 16
Description:Value stored to 'pedestal' during its initialization is never read

Annotated Source Code

1// $Id$
2//
3// File: DCDCHit_factory.cc
4// Created: Tue Aug 6 11:29:56 EDT 2013
5// Creator: davidl (on Darwin harriet.jlab.org 11.4.2 i386)
6//
7
8
9#include <iostream>
10#include <iomanip>
11using namespace std;
12
13#include <CDC/DCDCDigiHit.h>
14#include "DCDCHit_factory.h"
15#include "DCDCWire.h"
16#include <DAQ/Df125PulseIntegral.h>
17#include <DAQ/Df125Config.h>
18#include <DAQ/Df125CDCPulse.h>
19
20using namespace jana;
21
22static double DIGI_THRESHOLD = -1000000.0;
23//#define ENABLE_UPSAMPLING
24
25//------------------
26// init
27//------------------
28jerror_t DCDCHit_factory::init(void)
29{
30 gPARMS->SetDefaultParameter("CDC:DIGI_THRESHOLD",DIGI_THRESHOLD,
31 "Do not convert CDC digitized hits into DCDCHit objects"
32 " that would have q less than this");
33
34 // default values
35 Nrings = 0;
36 a_scale = 0.;
37 t_scale = 0.;
38 t_base = 0.;
39
40 // Set default number of number of detector channels
41 maxChannels = 3522;
42
43 /// set the base conversion scales
44 a_scale = 4.0E3/1.0E2;
45 t_scale = 8.0/10.0; // 8 ns/count and integer time is in 1/10th of sample
46 t_base = 0.; // ns
47
48 return NOERROR;
49}
50
51//------------------
52// brun
53//------------------
54jerror_t DCDCHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
55{
56 // Only print messages for one thread whenever run number change
57 static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
58 static set<int> runs_announced;
59 pthread_mutex_lock(&print_mutex);
60 bool print_messages = false;
61 if(runs_announced.find(runnumber) == runs_announced.end()){
62 print_messages = true;
63 runs_announced.insert(runnumber);
64 }
65 pthread_mutex_unlock(&print_mutex);
66
67 // calculate the number of straws in each ring
68 CalcNstraws(eventLoop, runnumber, Nstraws);
69 Nrings = Nstraws.size();
70
71 /// Read in calibration constants
72 vector<double> raw_gains;
73 vector<double> raw_pedestals;
74 vector<double> raw_time_offsets;
75
76 if(print_messages) jout << "In DCDCHit_factory, loading constants..." << std::endl;
77
78 // load scale factors
79 map<string,double> scale_factors;
80 if (eventLoop->GetCalib("/CDC/digi_scales", scale_factors))
81 jout << "Error loading /CDC/digi_scales !" << endl;
82 if (scale_factors.find("CDC_ADC_ASCALE") != scale_factors.end())
83 a_scale = scale_factors["CDC_ADC_ASCALE"];
84 else
85 jerr << "Unable to get CDC_ADC_ASCALE from /CDC/digi_scales !" << endl;
86
87#ifdef ENABLE_UPSAMPLING
88 //t_scale=1.;
89#else
90 if (scale_factors.find("CDC_ADC_TSCALE") != scale_factors.end())
91 t_scale = scale_factors["CDC_ADC_TSCALE"];
92 else
93 jerr << "Unable to get CDC_ADC_TSCALE from /CDC/digi_scales !" << endl;
94#endif
95
96 // load base time offset
97 map<string,double> base_time_offset;
98 if (eventLoop->GetCalib("/CDC/base_time_offset",base_time_offset))
99 jout << "Error loading /CDC/base_time_offset !" << endl;
100 if (base_time_offset.find("CDC_BASE_TIME_OFFSET") != base_time_offset.end())
101 t_base = base_time_offset["CDC_BASE_TIME_OFFSET"];
102 else
103 jerr << "Unable to get CDC_BASE_TIME_OFFSET from /CDC/base_time_offset !" << endl;
104
105 // load constant tables
106 if (eventLoop->GetCalib("/CDC/wire_gains", raw_gains))
107 jout << "Error loading /CDC/wire_gains !" << endl;
108 if (eventLoop->GetCalib("/CDC/pedestals", raw_pedestals))
109 jout << "Error loading /CDC/pedestals !" << endl;
110 if (eventLoop->GetCalib("/CDC/timing_offsets", raw_time_offsets))
111 jout << "Error loading /CDC/timing_offsets !" << endl;
112
113 // fill the tables
114 FillCalibTable(gains, raw_gains, Nstraws);
115 FillCalibTable(pedestals, raw_pedestals, Nstraws);
116 FillCalibTable(time_offsets, raw_time_offsets, Nstraws);
117
118 // Verify that the right number of rings was read for each set of constants
119 char str[256];
120 if (gains.size() != Nrings) {
121 sprintf(str, "Bad # of rings for CDC gain from CCDB! CCDB=%zu , should be %d", gains.size(), Nrings);
122 std::cerr << str << std::endl;
123 throw JException(str);
124 }
125 if (pedestals.size() != Nrings) {
126 sprintf(str, "Bad # of rings for CDC pedestal from CCDB! CCDB=%zu , should be %d", pedestals.size(), Nrings);
127 std::cerr << str << std::endl;
128 throw JException(str);
129 }
130 if (time_offsets.size() != Nrings) {
131 sprintf(str, "Bad # of rings for CDC time_offset from CCDB!"
132 " CCDB=%zu , should be %d", time_offsets.size(), Nrings);
133 std::cerr << str << std::endl;
134 throw JException(str);
135 }
136
137 // Verify the right number of straws was read for each ring for each set of constants
138 for (unsigned int i=0; i < Nrings; i++) {
139 if (gains[i].size() != Nstraws[i]) {
140 sprintf(str, "Bad # of straws for CDC gain from CCDB!"
141 " CCDB=%zu , should be %d for ring %d",
142 gains[i].size(), Nstraws[i], i+1);
143 std::cerr << str << std::endl;
144 throw JException(str);
145 }
146 if (pedestals[i].size() != Nstraws[i]) {
147 sprintf(str, "Bad # of straws for CDC pedestal from CCDB!"
148 " CCDB=%zu , should be %d for ring %d",
149 pedestals[i].size(), Nstraws[i], i+1);
150 std::cerr << str << std::endl;
151 throw JException(str);
152 }
153 if (time_offsets[i].size() != Nstraws[i]) {
154 sprintf(str, "Bad # of straws for CDC time_offset from CCDB!"
155 " CCDB=%zu , should be %d for ring %d",
156 time_offsets[i].size(), Nstraws[i], i+1);
157 std::cerr << str << std::endl;
158 throw JException(str);
159 }
160 }
161
162 return NOERROR;
163}
164
165//------------------
166// evnt
167//------------------
168jerror_t DCDCHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
169{
170 /// Generate DCDCHit object for each DCDCDigiHit object.
171 /// This is where the first set of calibration constants
172 /// is applied to convert from digitzed units into natural
173 /// units.
174 ///
175 /// Note that this code does NOT get called for simulated
176 /// data in HDDM format. The HDDM event source will copy
177 /// the precalibrated values directly into the _data vector.
178
179 /// In order to use the new Flash125 data types and maintain compatibility with the old code, what is below is a bit of a mess
180
181 vector<const DCDCDigiHit*> digihits;
182 loop->Get(digihits);
183 char str[256];
184 for (unsigned int i=0; i < digihits.size(); i++) {
185 const DCDCDigiHit *digihit = digihits[i];
186
187 if ( (digihit->QF & 0x1) != 0 ) continue; // Cut bad timing quality factor hits... (should check effect on efficiency)
188
189 const int &ring = digihit->ring;
190 const int &straw = digihit->straw;
191
192 // Make sure ring and straw are in valid range
193 if ( (ring < 1) || (ring > (int)Nrings)) {
194 sprintf(str, "DCDCDigiHit ring out of range!"
195 " ring=%d (should be 1-%d)", ring, Nrings);
196 throw JException(str);
197 }
198 if ( (straw < 1) || (straw > (int)Nstraws[ring-1])) {
199 sprintf(str, "DCDCDigiHit straw out of range!"
200 " straw=%d for ring=%d (should be 1-%d)",
201 straw, ring, Nstraws[ring-1]);
202 throw JException(str);
203 }
204
205 // Get pedestal. Preference is given to pedestal measured
206 // for event. Otherwise, use statistical one from CCDB
207 double pedestal = pedestals[ring-1][straw-1];
Value stored to 'pedestal' during its initialization is never read
208
209 // Grab the pedestal from the digihit since this should be consistent between the old and new formats
210 uint32_t raw_ped = digihit->pedestal;
211 uint32_t nsamples_integral = digihit->nsamples_integral;
212
213 // There are a few values from the new data type that are critical for the interpretation of the data
214 uint16_t IBIT = 0; // 2^{IBIT} Scale factor for integral
215 // uint16_t ABIT = 0; // 2^{ABIT} Scale factor for amplitude
216 uint16_t PBIT = 0; // 2^{PBIT} Scale factor for pedestal
217 uint16_t NW = 0;
218
219 // This is the place to make quality cuts on the data.
220 // Try to get the new data type, if that fails, try to get the old...
221 const Df125CDCPulse *CDCPulseObj = NULL__null;
222 digihit->GetSingle(CDCPulseObj);
223 if (CDCPulseObj != NULL__null){
224 const Df125Config *config = NULL__null;
225 CDCPulseObj->GetSingle(config);
226
227 // Set some constants to defaults until they appear correctly in the config words in the future
228 if(config){
229 IBIT = config->IBIT == 0xffff ? 4 : config->IBIT;
230 // ABIT = config->ABIT == 0xffff ? 3 : config->ABIT;
231 PBIT = config->PBIT == 0xffff ? 0 : config->PBIT;
232 NW = config->NW == 0xffff ? 180 : config->NW;
233 }else{
234 static int Nwarnings = 0;
235 if(Nwarnings<10){
236 _DBG_std::cerr<<"libraries/CDC/DCDCHit_factory.cc"<<":"
<<236<<" " << "NO Df125Config object associated with Df125FDCPulse object!" << endl;
237 Nwarnings++;
238 if(Nwarnings==10) _DBG_std::cerr<<"libraries/CDC/DCDCHit_factory.cc"<<":"
<<238<<" " << " --- LAST WARNING!! ---" << endl;
239 }
240 }
241 if(NW==0) NW=180; // some data was taken (<=run 4700) where NW was written as 0 to file
242
243 // The integration window in the CDC should always extend past the end of the window
244 // Only true after about run 4100
245 nsamples_integral = (NW - (digihit->pulse_time / 10));
246
247 }
248 else{ // Use the old format
249 // This code will at some point become deprecated in the future...
250 // This applies to the firmware for data taken until the fall of 2015.
251 // Mode 8: Raw data and processed data (except pulse integral).
252 // Mode 7: Processed data only.
253
254 // This error state is only present in mode 8
255 if (digihit->pulse_time==0.) continue;
256
257 // There is a slight difference between Mode 7 and 8 data
258 // The following condition signals an error state in the flash algorithm
259 // Do not make hits out of these
260 const Df125PulsePedestal* PPobj = NULL__null;
261 digihit->GetSingle(PPobj);
262 if (PPobj != NULL__null){
263 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
264 if (PPobj->pulse_number == 1) continue; // Unintentionally had 2 pulses found in fall 2014 data (0-1 counting issue)
265 }
266
267 const Df125PulseIntegral* PIobj = NULL__null;
268 digihit->GetSingle(PIobj);
269 if (PPobj == NULL__null || PIobj == NULL__null) continue; // We don't want hits where ANY of the associated information is missing
270 }
271
272 // Complete the pedestal subtracion here since we should know the correct number of samples.
273 uint32_t scaled_ped = raw_ped << PBIT;
274 pedestal = double(scaled_ped * nsamples_integral);
275
276 // Apply calibration constants here
277 double integral = double(digihit->pulse_integral << IBIT);
278 double t_raw = double(digihit->pulse_time);
279
280 double q = a_scale * gains[ring-1][straw-1] * (integral - pedestal);
281 double t = t_scale * t_raw - time_offsets[ring-1][straw-1] + t_base;
282
283 if (q < DIGI_THRESHOLD)
284 continue;
285
286 DCDCHit *hit = new DCDCHit;
287 hit->ring = ring;
288 hit->straw = straw;
289
290 // Values for d, itrack, ptype only apply to MC data
291 // note that ring/straw counting starts at 1
292 hit->q = q;
293 hit->t = t;
294 hit->d = 0.0;
295 hit->itrack = -1;
296 hit->ptype = 0;
297
298 hit->AddAssociatedObject(digihit);
299
300 _data.push_back(hit);
301 }
302
303 return NOERROR;
304}
305
306//------------------
307// erun
308//------------------
309jerror_t DCDCHit_factory::erun(void)
310{
311 return NOERROR;
312}
313
314//------------------
315// fini
316//------------------
317jerror_t DCDCHit_factory::fini(void)
318{
319 return NOERROR;
320}
321
322//------------------
323// CalcNstraws
324//------------------
325void DCDCHit_factory::CalcNstraws(jana::JEventLoop *eventLoop, int32_t runnumber, vector<unsigned int> &Nstraws)
326{
327 DGeometry *dgeom;
328 vector<vector<DCDCWire *> >cdcwires;
329
330 // Get pointer to DGeometry object
331 DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
332 dgeom = dapp->GetDGeometry(runnumber);
333
334 // Get the CDC wire table from the XML
335 dgeom->GetCDCWires(cdcwires);
336
337 // Fill array with the number of straws for each layer
338 // Also keep track of the total number of straws, i.e., the total number of detector channels
339 maxChannels = 0;
340 Nstraws.clear();
341 for (unsigned int i=0; i<cdcwires.size(); i++) {
342 Nstraws.push_back( cdcwires[i].size() );
343 maxChannels += cdcwires[i].size();
344 }
345
346 // clear up all of the wire information
347 for (unsigned int i=0; i<cdcwires.size(); i++) {
348 for (unsigned int j=0; j<cdcwires[i].size(); j++) {
349 delete cdcwires[i][j];
350 }
351 }
352 cdcwires.clear();
353}
354
355
356//------------------
357// FillCalibTable
358//------------------
359void DCDCHit_factory::FillCalibTable(vector< vector<double> > &table, vector<double> &raw_table,
360 vector<unsigned int> &Nstraws)
361{
362 int ring = 0;
363 int straw = 0;
364
365 // reset table before filling it
366 table.clear();
367 table.resize( Nstraws.size() );
368
369 for (unsigned int channel=0; channel<raw_table.size(); channel++,straw++) {
370 // make sure that we don't try to load info for channels that don't exist
371 if (channel == maxChannels) break;
372
373 // if we've hit the end of the ring, move on to the next
374 if (straw == (int)Nstraws[ring]) {
375 ring++;
376 straw = 0;
377 }
378
379 table[ring].push_back( raw_table[channel] );
380 }
381
382}
383
384
385//------------------------------------
386// GetConstant
387// Allow a few different interfaces
388//------------------------------------
389const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
390 const int in_ring, const int in_straw) const {
391
392 char str[256];
393
394 if ( (in_ring <= 0) || (static_cast<unsigned int>(in_ring) > Nrings)) {
395 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
396 " requested=%d , should be %ud", in_ring, Nrings);
397 std::cerr << str << std::endl;
398 throw JException(str);
399 }
400 if ( (in_straw <= 0) ||
401 (static_cast<unsigned int>(in_straw) > Nstraws[in_ring]))
402 {
403 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
404 " requested=%d , should be %ud", in_straw, Nstraws[in_ring]);
405 std::cerr << str << std::endl;
406 throw JException(str);
407 }
408
409 return the_table[in_ring-1][in_straw-1];
410}
411
412const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
413 const DCDCDigiHit *in_digihit) const {
414
415 char str[256];
416
417 if ( (in_digihit->ring <= 0) ||
418 (static_cast<unsigned int>(in_digihit->ring) > Nrings))
419 {
420 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
421 " requested=%d , should be %ud", in_digihit->ring, Nrings);
422 std::cerr << str << std::endl;
423 throw JException(str);
424 }
425 if ( (in_digihit->straw <= 0) ||
426 (static_cast<unsigned int>(in_digihit->straw) > Nstraws[in_digihit->ring]))
427 {
428 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
429 " requested=%d , should be %ud",
430 in_digihit->straw, Nstraws[in_digihit->ring]);
431 std::cerr << str << std::endl;
432 throw JException(str);
433 }
434
435 return the_table[in_digihit->ring-1][in_digihit->straw-1];
436}
437
438const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
439 const DCDCHit *in_hit) const {
440
441 char str[256];
442
443 if ( (in_hit->ring <= 0) || (static_cast<unsigned int>(in_hit->ring) > Nrings)) {
444 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
445 " requested=%d , should be %ud", in_hit->ring, Nrings);
446 std::cerr << str << std::endl;
447 throw JException(str);
448 }
449 if ( (in_hit->straw <= 0) ||
450 (static_cast<unsigned int>(in_hit->straw) > Nstraws[in_hit->ring])) {
451 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
452 " requested=%d , should be %ud",
453 in_hit->straw, Nstraws[in_hit->ring]);
454 std::cerr << str << std::endl;
455 throw JException(str);
456 }
457
458 return the_table[in_hit->ring-1][in_hit->straw-1];
459}
460/*
461 const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
462 const DTranslationTable *ttab,
463 const int in_rocid, const int in_slot, const int in_channel) const {
464
465 char str[256];
466
467 DTranslationTable::csc_t daq_index = { in_rocid, in_slot, in_channel };
468 DTranslationTable::DChannelInfo channel_info = ttab->GetDetectorIndex(daq_index);
469
470 if ( (channel_info.cdc.ring <= 0) || (channel_info.cdc.ring > Nrings)) {
471 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
472 " requested=%d , should be %ud", channel_info.cdc.ring, Nrings);
473 std::cerr << str << std::endl;
474 throw JException(str);
475 }
476 if ( (channel_info.cdc.straw <= 0) ||
477 (channel_info.cdc.straw > Nstraws[channel_info.cdc.ring]))
478 {
479 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
480 " requested=%d , should be %ud",
481 channel_info.cdc.ring, Nstraws[channel_info.cdc.straw]);
482 std::cerr << str << std::endl;
483 throw JException(str);
484 }
485
486 return the_table[channel_info.cdc.ring-1][channel_info.cdc.straw-1];
487 }
488 */