File: | libraries/CDC/DCDCHit_factory.cc |
Location: | line 207, 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 | 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 | //------------------ |
309 | jerror_t DCDCHit_factory::erun(void) |
310 | { |
311 | return NOERROR; |
312 | } |
313 | |
314 | //------------------ |
315 | // fini |
316 | //------------------ |
317 | jerror_t DCDCHit_factory::fini(void) |
318 | { |
319 | return NOERROR; |
320 | } |
321 | |
322 | //------------------ |
323 | // CalcNstraws |
324 | //------------------ |
325 | void 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 | //------------------ |
359 | void 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 | //------------------------------------ |
389 | const 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 | |
412 | const 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 | |
438 | const 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 | */ |