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