File: | libraries/CDC/DCDCHit_factory.cc |
Location: | line 205, column 16 |
Description: | Value stored to 'pedestal' during its initialization is never read |
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> |
11 | using 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 | |
20 | using namespace jana; |
21 | |
22 | static double DIGI_THRESHOLD = -1000000.0; |
23 | //#define ENABLE_UPSAMPLING |
24 | |
25 | //------------------ |
26 | // init |
27 | //------------------ |
28 | jerror_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 | //------------------ |
54 | jerror_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 | //------------------ |
168 | jerror_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 | //------------------ |
307 | jerror_t DCDCHit_factory::erun(void) |
308 | { |
309 | return NOERROR; |
310 | } |
311 | |
312 | //------------------ |
313 | // fini |
314 | //------------------ |
315 | jerror_t DCDCHit_factory::fini(void) |
316 | { |
317 | return NOERROR; |
318 | } |
319 | |
320 | //------------------ |
321 | // CalcNstraws |
322 | //------------------ |
323 | void 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 | //------------------ |
357 | void 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 | //------------------------------------ |
387 | const 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 | |
410 | const 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 | |
436 | const 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 | */ |