1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | #include <iostream> |
10 | #include <iomanip> |
11 | using namespace std; |
12 | |
13 | #include <FDC/DFDCGeometry.h> |
14 | #include <FDC/DFDCCathodeDigiHit.h> |
15 | #include <FDC/DFDCWireDigiHit.h> |
16 | #include "DFDCHit_factory.h" |
17 | #include <DAQ/Df125PulseIntegral.h> |
18 | #include <DAQ/Df125PulsePedestal.h> |
19 | #include <DAQ/Df125Config.h> |
20 | #include <DAQ/Df125FDCPulse.h> |
21 | using namespace jana; |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | jerror_t DFDCHit_factory::init(void) |
28 | { |
29 | |
30 | a_scale = 2.4E4/1.3E5; |
31 | t_scale = 8.0/10.0; |
32 | t_base = 0.; |
33 | fadc_t_base = 0.; |
34 | |
35 | return NOERROR; |
36 | } |
37 | |
38 | |
39 | |
40 | |
41 | jerror_t DFDCHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) |
42 | { |
43 | |
44 | static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
45 | static set<int> runs_announced; |
46 | pthread_mutex_lock(&print_mutex); |
47 | bool print_messages = false; |
48 | if(runs_announced.find(runnumber) == runs_announced.end()){ |
49 | print_messages = true; |
50 | runs_announced.insert(runnumber); |
51 | } |
52 | pthread_mutex_unlock(&print_mutex); |
53 | |
54 | |
55 | a_gains.clear(); |
56 | a_pedestals.clear(); |
57 | timing_offsets.clear(); |
58 | |
59 | |
60 | if(print_messages) jout << "In DFDCHit_factory, loading constants..." << endl; |
61 | |
62 | |
63 | map<string,double> scale_factors; |
64 | if(eventLoop->GetCalib("/FDC/digi_scales", scale_factors)) |
65 | jout << "Error loading /FDC/digi_scales !" << endl; |
66 | if( scale_factors.find("FDC_ADC_ASCALE") != scale_factors.end() ) { |
67 | a_scale = scale_factors["FDC_ADC_ASCALE"]; |
68 | } else { |
69 | jerr << "Unable to get FDC_ADC_ASCALE from /FDC/digi_scales !" << endl; |
70 | } |
71 | if( scale_factors.find("FDC_ADC_TSCALE") != scale_factors.end() ) { |
72 | t_scale = scale_factors["FDC_ADC_TSCALE"]; |
73 | } else { |
74 | jerr << "Unable to get FDC_ADC_TSCALE from /FDC/digi_scales !" << endl; |
75 | } |
76 | |
77 | |
78 | map<string,double> base_time_offset; |
79 | if (eventLoop->GetCalib("/FDC/base_time_offset",base_time_offset)) |
80 | jout << "Error loading /FDC/base_time_offset !" << endl; |
81 | if (base_time_offset.find("FDC_BASE_TIME_OFFSET") != base_time_offset.end()) |
82 | fadc_t_base = base_time_offset["FDC_BASE_TIME_OFFSET"]; |
83 | else |
84 | jerr << "Unable to get FDC_BASE_TIME_OFFSET from /FDC/base_time_offset !" << endl; |
85 | if (base_time_offset.find("FDC_TDC_BASE_TIME_OFFSET") != base_time_offset.end()) |
86 | t_base = base_time_offset["FDC_TDC_BASE_TIME_OFFSET"]; |
87 | else |
88 | jerr << "Unable to get FDC_TDC_BASE_TIME_OFFSET from /FDC/base_time_offset !" << endl; |
89 | |
90 | |
91 | |
92 | LoadPackageCalibTables(eventLoop,"/FDC/package1"); |
93 | LoadPackageCalibTables(eventLoop,"/FDC/package2"); |
94 | LoadPackageCalibTables(eventLoop,"/FDC/package3"); |
95 | LoadPackageCalibTables(eventLoop,"/FDC/package4"); |
96 | |
97 | |
98 | char str[256]; |
99 | if(a_gains.size() != FDC_NUM_PLANES72) { |
100 | sprintf(str, "Bad # of planes for FDC gains from CCDB! CCDB=%zu , should be %d", |
101 | a_gains.size(), FDC_NUM_PLANES72); |
102 | cerr << str << endl; |
103 | throw JException(str); |
104 | } |
105 | if(a_pedestals.size() != FDC_NUM_PLANES72) { |
106 | sprintf(str, "Bad # of planes for FDC pedestals from CCDB! CCDB=%zu , should be %d", |
107 | a_pedestals.size(), FDC_NUM_PLANES72); |
108 | cerr << str << endl; |
109 | throw JException(str); |
110 | } |
111 | if(timing_offsets.size() != FDC_NUM_PLANES72) { |
112 | sprintf(str, "Bad # of planes for FDC timing offsets from CCDB! CCDB=%zu , should be %d", |
113 | timing_offsets.size(), FDC_NUM_PLANES72); |
114 | cerr << str << endl; |
115 | throw JException(str); |
116 | } |
117 | |
118 | return NOERROR; |
119 | } |
120 | |
121 | |
122 | |
123 | |
124 | jerror_t DFDCHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) |
125 | { |
126 | |
127 | |
128 | |
129 | |
130 | |
131 | |
132 | |
133 | |
134 | |
135 | char str[256]; |
136 | |
137 | const DTTabUtilities* locTTabUtilities = NULL__null; |
138 | loop->GetSingle(locTTabUtilities); |
139 | |
140 | |
141 | vector<const DFDCCathodeDigiHit*> cathodedigihits; |
142 | loop->Get(cathodedigihits); |
143 | for(unsigned int i=0; i<cathodedigihits.size(); i++){ |
144 | const DFDCCathodeDigiHit *digihit = cathodedigihits[i]; |
145 | |
146 | |
147 | |
148 | |
149 | |
150 | |
151 | |
152 | |
153 | |
154 | |
155 | |
156 | |
157 | |
158 | |
159 | |
160 | |
161 | |
162 | |
163 | |
164 | int layer=digihit->view; |
165 | int gLayer=digihit->chamber + 6*(digihit->package - 1); |
166 | int gPlane=layer + 3*(gLayer - 1); |
167 | |
168 | if((gPlane < 1) || (gPlane > FDC_NUM_PLANES72)) { |
169 | sprintf(str, "DFDCDigiHit plane out of range! gPlane=%d (should be 1-%d)", gPlane, FDC_NUM_PLANES72); |
170 | throw JException(str); |
171 | } |
172 | if( (digihit->strip < 1) || (digihit->strip > STRIPS_PER_PLANE192)) { |
173 | sprintf(str, "DFDCDigiHit straw out of range! strip=%d for plane=%d (should be 1-%d)", digihit->strip, gPlane, STRIPS_PER_PLANE192); |
174 | throw JException(str); |
175 | } |
176 | |
177 | int plane_index=gPlane-1; |
178 | int strip_index=digihit->strip-1; |
179 | |
180 | |
181 | double T = (double)digihit->pulse_time; |
182 | |
183 | |
184 | |
185 | double pedestal = a_pedestals[plane_index][strip_index]; |
| Value stored to 'pedestal' during its initialization is never read |
186 | |
187 | |
188 | uint32_t raw_ped = digihit->pedestal; |
189 | uint32_t nsamples_integral = digihit->nsamples_integral; |
190 | |
191 | |
192 | uint16_t IBIT = 0; |
193 | uint16_t ABIT = 0; |
194 | uint16_t PBIT = 0; |
195 | |
196 | |
197 | |
198 | uint32_t pulse_peak = 0; |
199 | const Df125FDCPulse *FDCPulseObj = NULL__null; |
200 | digihit->GetSingle(FDCPulseObj); |
201 | if (FDCPulseObj != NULL__null){ |
202 | |
203 | const Df125Config *config = NULL__null; |
204 | FDCPulseObj->GetSingle(config); |
205 | |
206 | |
207 | |
208 | IBIT = config->IBIT == 0xffff ? 4 : config->IBIT; |
209 | ABIT = config->ABIT == 0xffff ? 3 : config->ABIT; |
210 | PBIT = config->PBIT == 0xffff ? 0 : config->PBIT; |
211 | uint16_t NW = config->NW == 0xffff ? 80 : config->NW; |
212 | uint16_t IE = config->IE == 0xffff ? 16 : config->IE; |
213 | |
214 | if ((NW - (digihit->pulse_time / 10)) < IE){ |
215 | nsamples_integral = (NW - (digihit->pulse_time / 10)); |
216 | } |
217 | else{ |
218 | nsamples_integral = IE; |
219 | } |
220 | |
221 | pulse_peak = FDCPulseObj->peak_amp << ABIT; |
222 | } |
223 | else{ |
224 | |
225 | |
226 | |
227 | const Df125PulsePedestal* PPobj = NULL__null; |
228 | digihit->GetSingle(PPobj); |
229 | if (PPobj != NULL__null){ |
230 | if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue; |
231 | if (PPobj->pulse_number == 1) continue; |
232 | pulse_peak = PPobj->pulse_peak; |
233 | } |
234 | |
235 | const Df125PulseIntegral* PIobj = NULL__null; |
236 | digihit->GetSingle(PIobj); |
237 | if ( PPobj == NULL__null || PIobj == NULL__null) continue; |
238 | } |
239 | |
240 | |
241 | uint32_t scaled_ped = raw_ped << PBIT; |
242 | pedestal = double(scaled_ped * nsamples_integral); |
243 | |
244 | double integral = double(digihit->pulse_integral << IBIT); |
245 | |
246 | |
247 | |
248 | double q = a_scale * a_gains[plane_index][strip_index] * (integral-pedestal); |
249 | double t = t_scale * T - timing_offsets[plane_index][strip_index]+fadc_t_base; |
250 | |
251 | DFDCHit *hit = new DFDCHit; |
252 | hit->layer = digihit->view; |
253 | hit->gLayer = gLayer; |
254 | hit->gPlane = gPlane; |
255 | hit->module = 1 + (gLayer-1)/3; |
256 | hit->element = digihit->strip; |
257 | hit->plane = digihit->view; |
258 | hit->r = DFDCGeometry::getWireR(hit); |
259 | hit->d = 0.0; |
260 | hit->type = digihit->strip_type; |
261 | hit->itrack = -1; |
262 | hit->ptype = 0; |
263 | hit->q = q; |
264 | hit->t = t; |
265 | hit->pulse_height=a_gains[plane_index][strip_index] |
266 | *double(pulse_peak - scaled_ped); |
267 | |
268 | |
269 | |
270 | hit->AddAssociatedObject(digihit); |
271 | |
272 | _data.push_back(hit); |
273 | } |
274 | |
275 | |
276 | vector<const DFDCWireDigiHit*> wiredigihits; |
277 | loop->Get(wiredigihits); |
278 | for(unsigned int i=0; i<wiredigihits.size(); i++){ |
279 | const DFDCWireDigiHit *digihit = wiredigihits[i]; |
280 | |
281 | |
282 | |
283 | |
284 | |
285 | |
286 | |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | DFDCHit *hit = new DFDCHit; |
300 | hit->gLayer = digihit->chamber + 6*(digihit->package - 1); |
301 | hit->module = 1 + (hit->gLayer-1)/3; |
302 | hit->layer = hit->gLayer - (hit->module-1)*3; |
303 | hit->plane = 2; |
304 | hit->gPlane = hit->plane + 3*(hit->gLayer - 1); |
305 | hit->element = digihit->wire; |
306 | hit->r = DFDCGeometry::getWireR(hit); |
307 | hit->d = 0.0; |
308 | hit->type = DFDCHit::AnodeWire; |
309 | hit->itrack = -1; |
310 | hit->ptype = 0; |
311 | |
312 | |
313 | if( (hit->gPlane < 1) || (hit->gPlane > FDC_NUM_PLANES72)) { |
314 | sprintf(str, "DFDCDigiHit plane out of range! gPlane=%d (should be 1-%d)", hit->gPlane, FDC_NUM_PLANES72); |
315 | throw JException(str); |
316 | } |
317 | if( (digihit->wire < 1) || (digihit->wire > WIRES_PER_PLANE96)) { |
318 | sprintf(str, "DFDCDigiHit straw out of range! wire=%d for plane=%d (should be 1-%d)", digihit->wire, hit->gPlane, WIRES_PER_PLANE96); |
319 | throw JException(str); |
320 | } |
321 | |
322 | |
323 | double T = locTTabUtilities->Convert_DigiTimeToNs_F1TDC(digihit) - timing_offsets[hit->gPlane-1][hit->element-1] + t_base; |
324 | hit->q = 0.0; |
325 | hit->pulse_height=0.0; |
326 | hit->t = T; |
327 | |
328 | hit->AddAssociatedObject(digihit); |
329 | |
330 | _data.push_back(hit); |
331 | } |
332 | |
333 | return NOERROR; |
334 | } |
335 | |
336 | |
337 | |
338 | |
339 | jerror_t DFDCHit_factory::erun(void) |
340 | { |
341 | return NOERROR; |
342 | } |
343 | |
344 | |
345 | |
346 | |
347 | jerror_t DFDCHit_factory::fini(void) |
348 | { |
349 | return NOERROR; |
350 | } |
351 | |
352 | |
353 | |
354 | |
355 | |
356 | void DFDCHit_factory::LoadPackageCalibTables(jana::JEventLoop *eventLoop, string ccdb_prefix) |
357 | { |
358 | vector< vector<double> > new_gains, new_pedestals, new_strip_t0s, new_wire_t0s; |
359 | char str[256]; |
360 | |
361 | if(eventLoop->GetCalib(ccdb_prefix+"/strip_gains", new_gains)) |
362 | cout << "Error loading "+ccdb_prefix+"/strip_gains !" << endl; |
363 | if(eventLoop->GetCalib(ccdb_prefix+"/strip_pedestals", new_pedestals)) |
364 | cout << "Error loading "+ccdb_prefix+"/strip_pedestals !" << endl; |
365 | if(eventLoop->GetCalib(ccdb_prefix+"/strip_timing_offsets", new_strip_t0s)) |
366 | cout << "Error loading "+ccdb_prefix+"/strip_timing_offsets!" << endl; |
367 | if(eventLoop->GetCalib(ccdb_prefix+"/wire_timing_offsets", new_wire_t0s)) |
368 | cout << "Error loading "+ccdb_prefix+"/wire_timing_offsets!" << endl; |
369 | |
370 | for(int nchamber=0; nchamber<6; nchamber++) { |
371 | |
372 | |
373 | if(new_gains[2*nchamber].size() != STRIPS_PER_PLANE192) { |
374 | sprintf(str, "Bad # of strips for FDC gain from CCDB! CCDB=%zu , should be %d", new_gains[2*nchamber].size(), STRIPS_PER_PLANE192); |
375 | cerr << str << endl; |
376 | throw JException(str); |
377 | } |
378 | if(new_gains[2*nchamber+1].size() != STRIPS_PER_PLANE192) { |
379 | sprintf(str, "Bad # of strips for FDC gain from CCDB! CCDB=%zu , should be %d", new_gains[2*nchamber+1].size(), STRIPS_PER_PLANE192); |
380 | cerr << str << endl; |
381 | throw JException(str); |
382 | } |
383 | if(new_pedestals[2*nchamber].size() != STRIPS_PER_PLANE192) { |
384 | sprintf(str, "Bad # of strips for FDC pedestals from CCDB! CCDB=%zu , should be %d", new_pedestals[2*nchamber].size(), STRIPS_PER_PLANE192); |
385 | cerr << str << endl; |
386 | throw JException(str); |
387 | } |
388 | if(new_pedestals[2*nchamber+1].size() != STRIPS_PER_PLANE192) { |
389 | sprintf(str, "Bad # of strips for FDC pedestals from CCDB! CCDB=%zu , should be %d", new_pedestals[2*nchamber+1].size(), STRIPS_PER_PLANE192); |
390 | cerr << str << endl; |
391 | throw JException(str); |
392 | } |
393 | if(new_strip_t0s[2*nchamber].size() != STRIPS_PER_PLANE192) { |
394 | sprintf(str, "Bad # of strips for FDC timing offsets from CCDB! CCDB=%zu , should be %d", new_strip_t0s[2*nchamber].size(), STRIPS_PER_PLANE192); |
395 | cerr << str << endl; |
396 | throw JException(str); |
397 | } |
398 | if(new_strip_t0s[2*nchamber+1].size() != STRIPS_PER_PLANE192) { |
399 | sprintf(str, "Bad # of strips for FDC timing offsets from CCDB! CCDB=%zu , should be %d", new_strip_t0s[2*nchamber+1].size(), STRIPS_PER_PLANE192); |
400 | cerr << str << endl; |
401 | throw JException(str); |
402 | } |
403 | if(new_wire_t0s[nchamber].size() != WIRES_PER_PLANE96) { |
404 | sprintf(str, "Bad # of wires for FDC timing offsets from CCDB! CCDB=%zu , should be %d", new_wire_t0s[2*nchamber].size(), WIRES_PER_PLANE96); |
405 | cerr << str << endl; |
406 | throw JException(str); |
407 | } |
408 | |
409 | |
410 | |
411 | a_gains.push_back( new_gains[2*nchamber] ); |
412 | a_gains.push_back( vector<double>() ); |
413 | a_gains.push_back( new_gains[2*nchamber+1] ); |
414 | |
415 | |
416 | a_pedestals.push_back( new_pedestals[2*nchamber] ); |
417 | a_pedestals.push_back( vector<double>() ); |
418 | a_pedestals.push_back( new_pedestals[2*nchamber+1] ); |
419 | |
420 | |
421 | timing_offsets.push_back( new_strip_t0s[2*nchamber] ); |
422 | timing_offsets.push_back( new_wire_t0s[nchamber] ); |
423 | timing_offsets.push_back( new_strip_t0s[2*nchamber+1] ); |
424 | |
425 | } |
426 | } |
427 | |
428 | |
429 | |
430 | |
431 | |
432 | const double DFDCHit_factory::GetConstant(const fdc_digi_constants_t &the_table, |
433 | const int in_gPlane, const int in_element) const { |
434 | |
435 | char str[256]; |
436 | |
437 | if( (in_gPlane <= 0) || (static_cast<unsigned int>(in_gPlane) > FDC_NUM_PLANES72)) { |
438 | sprintf(str, "Bad gPlane # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %ud", in_gPlane, FDC_NUM_PLANES72); |
439 | cerr << str << endl; |
440 | throw JException(str); |
441 | } |
442 | |
443 | if( (in_element <= 0) || (static_cast<unsigned int>(in_element) > the_table[in_gPlane].size())) { |
444 | sprintf(str, "Bad element # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %zu", in_element, the_table[in_gPlane].size()); |
445 | cerr << str << endl; |
446 | throw JException(str); |
447 | } |
448 | |
449 | return the_table[in_gPlane-1][in_element-1]; |
450 | } |
451 | |
452 | const double DFDCHit_factory::GetConstant(const fdc_digi_constants_t &the_table, |
453 | const DFDCCathodeDigiHit *in_digihit) const { |
454 | |
455 | char str[256]; |
456 | |
457 | int gLayer = in_digihit->chamber + 6*(in_digihit->package - 1); |
458 | int gPlane = in_digihit->view + 3*(gLayer - 1); |
459 | |
460 | if( (gPlane <= 0) || (static_cast<unsigned int>(gPlane) > FDC_NUM_PLANES72)) { |
461 | sprintf(str, "Bad gPlane # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %ud", gPlane, FDC_NUM_PLANES72); |
462 | cerr << str << endl; |
463 | throw JException(str); |
464 | } |
465 | |
466 | if( (in_digihit->strip <= 0) || (static_cast<unsigned int>(in_digihit->strip) > STRIPS_PER_PLANE192)) { |
467 | sprintf(str, "Bad strip # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %ud", in_digihit->strip, STRIPS_PER_PLANE192); |
468 | cerr << str << endl; |
469 | throw JException(str); |
470 | } |
471 | |
472 | return the_table[gPlane-1][in_digihit->strip-1]; |
473 | } |
474 | |
475 | const double DFDCHit_factory::GetConstant(const fdc_digi_constants_t &the_table, |
476 | const DFDCWireDigiHit *in_digihit) const { |
477 | |
478 | char str[256]; |
479 | |
480 | int gLayer = in_digihit->chamber + 6*(in_digihit->package - 1); |
481 | int gPlane = 2 + 3*(gLayer - 1); |
482 | |
483 | if( (gPlane <= 0) || (static_cast<unsigned int>(gPlane) > FDC_NUM_PLANES72)) { |
484 | sprintf(str, "Bad gPlane # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %ud", gPlane, FDC_NUM_PLANES72); |
485 | cerr << str << endl; |
486 | throw JException(str); |
487 | } |
488 | |
489 | if( (in_digihit->wire <= 0) || (static_cast<unsigned int>(in_digihit->wire) > WIRES_PER_PLANE96)) { |
490 | sprintf(str, "Bad wire # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %ud", in_digihit->wire, WIRES_PER_PLANE96); |
491 | cerr << str << endl; |
492 | throw JException(str); |
493 | } |
494 | |
495 | return the_table[gPlane-1][in_digihit->wire-1]; |
496 | } |
497 | |
498 | const double DFDCHit_factory::GetConstant(const fdc_digi_constants_t &the_table, |
499 | const DFDCHit *in_hit) const { |
500 | |
501 | char str[256]; |
502 | |
503 | if( (in_hit->gPlane <= 0) || (static_cast<unsigned int>(in_hit->gPlane) > FDC_NUM_PLANES72)) { |
504 | sprintf(str, "Bad gPlane # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %ud", in_hit->gPlane, FDC_NUM_PLANES72); |
505 | cerr << str << endl; |
506 | throw JException(str); |
507 | } |
508 | |
509 | if( (in_hit->element <= 0) || (static_cast<unsigned int>(in_hit->element) > the_table[in_hit->gPlane].size())) { |
510 | sprintf(str, "Bad element # requested in DFDCHit_factory::GetConstant()! requested=%d , should be %zu", in_hit->element, the_table[in_hit->gPlane].size()); |
511 | cerr << str << endl; |
512 | throw JException(str); |
513 | } |
514 | |
515 | return the_table[in_hit->gPlane-1][in_hit->element-1]; |
516 | } |
517 | |
518 | |
519 | |
520 | |
521 | |
522 | |
523 | |
524 | |
525 | |
526 | |
527 | |
528 | |
529 | |
530 | |
531 | |
532 | |
533 | |
534 | |
535 | |
536 | |
537 | |
538 | |
539 | |
540 | |
541 | |
542 | |
543 | |
544 | |
545 | |
546 | |
547 | |
548 | |
549 | |
550 | |
551 | |
552 | |
553 | |
554 | |
555 | |
556 | |
557 | |
558 | |
559 | |
560 | |
561 | |
562 | |
563 | |
564 | |
565 | |
566 | |
567 | |
568 | |
569 | |
570 | |
571 | |
572 | |
573 | |