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