Bug Summary

File:libraries/CDC/DCDCHit_factory.cc
Location:line 205, 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 const int &ring = digihit->ring;
188 const int &straw = digihit->straw;
189
190 // Make sure ring and straw are in valid range
191 if ( (ring < 1) || (ring > (int)Nrings)) {
192 sprintf(str, "DCDCDigiHit ring out of range!"
193 " ring=%d (should be 1-%d)", ring, Nrings);
194 throw JException(str);
195 }
196 if ( (straw < 1) || (straw > (int)Nstraws[ring-1])) {
197 sprintf(str, "DCDCDigiHit straw out of range!"
198 " straw=%d for ring=%d (should be 1-%d)",
199 straw, ring, Nstraws[ring-1]);
200 throw JException(str);
201 }
202
203 // Get pedestal. Preference is given to pedestal measured
204 // for event. Otherwise, use statistical one from CCDB
205 double pedestal = pedestals[ring-1][straw-1];
Value stored to 'pedestal' during its initialization is never read
206
207 // Grab the pedestal from the digihit since this should be consistent between the old and new formats
208 uint32_t raw_ped = digihit->pedestal;
209 uint32_t nsamples_integral = digihit->nsamples_integral;
210
211 // There are a few values from the new data type that are critical for the interpretation of the data
212 uint16_t IBIT = 0; // 2^{IBIT} Scale factor for integral
213 uint16_t ABIT = 0; // 2^{ABIT} Scale factor for amplitude
214 uint16_t PBIT = 0; // 2^{PBIT} Scale factor for pedestal
215
216 // This is the place to make quality cuts on the data.
217 // Try to get the new data type, if that fails, try to get the old...
218 const Df125CDCPulse *CDCPulseObj = NULL__null;
219 digihit->GetSingle(CDCPulseObj);
220 if (CDCPulseObj != NULL__null){
221 const Df125Config *config = NULL__null;
222 CDCPulseObj->GetSingle(config);
223
224 // Set some constants to defaults until they appear correctly in the config words in the future
225 IBIT = config->IBIT == 0xffff ? 4 : config->IBIT;
226 ABIT = config->ABIT == 0xffff ? 3 : config->ABIT;
227 PBIT = config->PBIT == 0xffff ? 0 : config->PBIT;
228 uint16_t NW = config->NW == 0xffff ? 200 : config->NW;
229
230 // The integration window in the CDC should always extend past the end of the window
231 // Only true after about run 4100
232 nsamples_integral = (NW - (digihit->pulse_time / 10));
233
234 }
235 else{ // Use the old format
236 // This code will at some point become deprecated in the future...
237 // This applies to the firmware for data taken until the fall of 2015.
238 // Mode 8: Raw data and processed data (except pulse integral).
239 // Mode 7: Processed data only.
240
241 // This error state is only present in mode 8
242 if (digihit->pulse_time==0.) continue;
243
244 // There is a slight difference between Mode 7 and 8 data
245 // The following condition signals an error state in the flash algorithm
246 // Do not make hits out of these
247 const Df125PulsePedestal* PPobj = NULL__null;
248 digihit->GetSingle(PPobj);
249 if (PPobj != NULL__null){
250 if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue;
251 if (PPobj->pulse_number == 1) continue; // Unintentionally had 2 pulses found in fall 2014 data (0-1 counting issue)
252 }
253
254 const Df125PulseIntegral* PIobj = NULL__null;
255 digihit->GetSingle(PIobj);
256 if (PPobj == NULL__null || PIobj == NULL__null) continue; // We don't want hits where ANY of the associated information is missing
257 }
258
259 // Complete the pedestal subtracion here since we should know the correct number of samples.
260 uint32_t scaled_ped = raw_ped << PBIT;
261 pedestal = double(scaled_ped * nsamples_integral);
262
263 // Apply calibration constants here
264 double integral = double(digihit->pulse_integral << IBIT);
265 double t_raw = double(digihit->pulse_time);
266
267 double q = a_scale * gains[ring-1][straw-1] * (integral - pedestal);
268 double t = t_scale * t_raw - time_offsets[ring-1][straw-1] + t_base;
269
270 if (q < DIGI_THRESHOLD)
271 continue;
272
273 DCDCHit *hit = new DCDCHit;
274 hit->ring = ring;
275 hit->straw = straw;
276
277 // Values for d, itrack, ptype only apply to MC data
278 // note that ring/straw counting starts at 1
279 hit->q = q;
280 hit->t = t;
281 hit->d = 0.0;
282 hit->itrack = -1;
283 hit->ptype = 0;
284
285 hit->AddAssociatedObject(digihit);
286
287 _data.push_back(hit);
288 }
289
290 return NOERROR;
291}
292
293//------------------
294// erun
295//------------------
296jerror_t DCDCHit_factory::erun(void)
297{
298 return NOERROR;
299}
300
301//------------------
302// fini
303//------------------
304jerror_t DCDCHit_factory::fini(void)
305{
306 return NOERROR;
307}
308
309//------------------
310// CalcNstraws
311//------------------
312void DCDCHit_factory::CalcNstraws(jana::JEventLoop *eventLoop, int32_t runnumber, vector<unsigned int> &Nstraws)
313{
314 DGeometry *dgeom;
315 vector<vector<DCDCWire *> >cdcwires;
316
317 // Get pointer to DGeometry object
318 DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
319 dgeom = dapp->GetDGeometry(runnumber);
320
321 // Get the CDC wire table from the XML
322 dgeom->GetCDCWires(cdcwires);
323
324 // Fill array with the number of straws for each layer
325 // Also keep track of the total number of straws, i.e., the total number of detector channels
326 maxChannels = 0;
327 Nstraws.clear();
328 for (unsigned int i=0; i<cdcwires.size(); i++) {
329 Nstraws.push_back( cdcwires[i].size() );
330 maxChannels += cdcwires[i].size();
331 }
332
333 // clear up all of the wire information
334 for (unsigned int i=0; i<cdcwires.size(); i++) {
335 for (unsigned int j=0; j<cdcwires[i].size(); j++) {
336 delete cdcwires[i][j];
337 }
338 }
339 cdcwires.clear();
340}
341
342
343//------------------
344// FillCalibTable
345//------------------
346void DCDCHit_factory::FillCalibTable(vector< vector<double> > &table, vector<double> &raw_table,
347 vector<unsigned int> &Nstraws)
348{
349 int ring = 0;
350 int straw = 0;
351
352 // reset table before filling it
353 table.clear();
354 table.resize( Nstraws.size() );
355
356 for (unsigned int channel=0; channel<raw_table.size(); channel++,straw++) {
357 // make sure that we don't try to load info for channels that don't exist
358 if (channel == maxChannels) break;
359
360 // if we've hit the end of the ring, move on to the next
361 if (straw == (int)Nstraws[ring]) {
362 ring++;
363 straw = 0;
364 }
365
366 table[ring].push_back( raw_table[channel] );
367 }
368
369}
370
371
372//------------------------------------
373// GetConstant
374// Allow a few different interfaces
375//------------------------------------
376const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
377 const int in_ring, const int in_straw) const {
378
379 char str[256];
380
381 if ( (in_ring <= 0) || (static_cast<unsigned int>(in_ring) > Nrings)) {
382 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
383 " requested=%d , should be %ud", in_ring, Nrings);
384 std::cerr << str << std::endl;
385 throw JException(str);
386 }
387 if ( (in_straw <= 0) ||
388 (static_cast<unsigned int>(in_straw) > Nstraws[in_ring]))
389 {
390 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
391 " requested=%d , should be %ud", in_straw, Nstraws[in_ring]);
392 std::cerr << str << std::endl;
393 throw JException(str);
394 }
395
396 return the_table[in_ring-1][in_straw-1];
397}
398
399const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
400 const DCDCDigiHit *in_digihit) const {
401
402 char str[256];
403
404 if ( (in_digihit->ring <= 0) ||
405 (static_cast<unsigned int>(in_digihit->ring) > Nrings))
406 {
407 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
408 " requested=%d , should be %ud", in_digihit->ring, Nrings);
409 std::cerr << str << std::endl;
410 throw JException(str);
411 }
412 if ( (in_digihit->straw <= 0) ||
413 (static_cast<unsigned int>(in_digihit->straw) > Nstraws[in_digihit->ring]))
414 {
415 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
416 " requested=%d , should be %ud",
417 in_digihit->straw, Nstraws[in_digihit->ring]);
418 std::cerr << str << std::endl;
419 throw JException(str);
420 }
421
422 return the_table[in_digihit->ring-1][in_digihit->straw-1];
423}
424
425const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
426 const DCDCHit *in_hit) const {
427
428 char str[256];
429
430 if ( (in_hit->ring <= 0) || (static_cast<unsigned int>(in_hit->ring) > Nrings)) {
431 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
432 " requested=%d , should be %ud", in_hit->ring, Nrings);
433 std::cerr << str << std::endl;
434 throw JException(str);
435 }
436 if ( (in_hit->straw <= 0) ||
437 (static_cast<unsigned int>(in_hit->straw) > Nstraws[in_hit->ring])) {
438 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
439 " requested=%d , should be %ud",
440 in_hit->straw, Nstraws[in_hit->ring]);
441 std::cerr << str << std::endl;
442 throw JException(str);
443 }
444
445 return the_table[in_hit->ring-1][in_hit->straw-1];
446}
447/*
448 const double DCDCHit_factory::GetConstant(const cdc_digi_constants_t &the_table,
449 const DTranslationTable *ttab,
450 const int in_rocid, const int in_slot, const int in_channel) const {
451
452 char str[256];
453
454 DTranslationTable::csc_t daq_index = { in_rocid, in_slot, in_channel };
455 DTranslationTable::DChannelInfo channel_info = ttab->GetDetectorIndex(daq_index);
456
457 if ( (channel_info.cdc.ring <= 0) || (channel_info.cdc.ring > Nrings)) {
458 sprintf(str, "Bad ring # requested in DCDCHit_factory::GetConstant()!"
459 " requested=%d , should be %ud", channel_info.cdc.ring, Nrings);
460 std::cerr << str << std::endl;
461 throw JException(str);
462 }
463 if ( (channel_info.cdc.straw <= 0) ||
464 (channel_info.cdc.straw > Nstraws[channel_info.cdc.ring]))
465 {
466 sprintf(str, "Bad straw # requested in DCDCHit_factory::GetConstant()!"
467 " requested=%d , should be %ud",
468 channel_info.cdc.ring, Nstraws[channel_info.cdc.straw]);
469 std::cerr << str << std::endl;
470 throw JException(str);
471 }
472
473 return the_table[channel_info.cdc.ring-1][channel_info.cdc.straw-1];
474 }
475 */