1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | #include <iostream> |
10 | #include <iomanip> |
11 | #include <cmath> |
12 | #include <vector> |
13 | #include <limits> |
14 | |
15 | using namespace std; |
16 | |
17 | #include <TOF/DTOFDigiHit.h> |
18 | #include <TOF/DTOFTDCDigiHit.h> |
19 | #include "DTOFHit_factory.h" |
20 | #include <DAQ/Df250PulseIntegral.h> |
21 | #include <DAQ/Df250Config.h> |
22 | #include <DAQ/DCODAROCInfo.h> |
23 | using namespace jana; |
24 | |
25 | static bool COSMIC_DATA = false; |
26 | |
27 | |
28 | |
29 | |
30 | jerror_t DTOFHit_factory::init(void) |
31 | { |
32 | DELTA_T_ADC_TDC_MAX = 10.0; |
33 | |
34 | gPARMS->SetDefaultParameter("TOF:DELTA_T_ADC_TDC_MAX", DELTA_T_ADC_TDC_MAX, |
35 | "Maximum difference in ns between a (calibrated) fADC time and F1TDC time for them to be matched in a single hit"); |
36 | |
37 | int analyze_cosmic_data = 0; |
38 | gPARMS->SetDefaultParameter("TOF:COSMIC_DATA", analyze_cosmic_data, |
39 | "Special settings for analysing cosmic data"); |
40 | if(analyze_cosmic_data > 0) |
41 | COSMIC_DATA = true; |
42 | |
43 | CHECK_FADC_ERRORS = false; |
44 | gPARMS->SetDefaultParameter("TOF:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits"); |
45 | |
46 | |
47 | |
48 | a_scale = 0.2/5.2E5; |
49 | t_scale = 0.0625; |
50 | t_base = 0.; |
51 | t_base_tdc = 0.; |
52 | |
53 | if(COSMIC_DATA) |
54 | |
55 | tdc_adc_time_offset = 110.; |
56 | else |
57 | tdc_adc_time_offset = 0.; |
58 | |
59 | TOF_NUM_PLANES = 2; |
60 | TOF_NUM_BARS = 44; |
61 | |
62 | return NOERROR; |
63 | } |
64 | |
65 | |
66 | |
67 | |
68 | jerror_t DTOFHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) |
69 | { |
70 | |
71 | static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
72 | static set<int> runs_announced; |
73 | pthread_mutex_lock(&print_mutex); |
74 | bool print_messages = false; |
75 | if(runs_announced.find(runnumber) == runs_announced.end()){ |
76 | print_messages = true; |
77 | runs_announced.insert(runnumber); |
78 | } |
79 | pthread_mutex_unlock(&print_mutex); |
80 | |
81 | |
82 | vector<const DTOFGeometry*> tofGeomVect; |
83 | eventLoop->Get( tofGeomVect ); |
84 | if(tofGeomVect.size()<1) return OBJECT_NOT_AVAILABLE; |
85 | const DTOFGeometry& tofGeom = *(tofGeomVect[0]); |
86 | |
87 | |
88 | vector<double> raw_adc_pedestals; |
89 | vector<double> raw_adc_gains; |
90 | vector<double> raw_adc_offsets; |
91 | vector<double> raw_tdc_offsets; |
92 | vector<double> raw_adc2E; |
93 | |
94 | if(print_messages) jout << "In DTOFHit_factory, loading constants..." << endl; |
95 | |
96 | |
97 | map<string,double> scale_factors; |
98 | if(eventLoop->GetCalib("/TOF/digi_scales", scale_factors)) |
99 | jout << "Error loading /TOF/digi_scales !" << endl; |
100 | if( scale_factors.find("TOF_ADC_ASCALE") != scale_factors.end() ) { |
101 | ; |
102 | } else { |
103 | jerr << "Unable to get TOF_ADC_ASCALE from /TOF/digi_scales !" << endl; |
104 | } |
105 | if( scale_factors.find("TOF_ADC_TSCALE") != scale_factors.end() ) { |
106 | ; |
107 | } else { |
108 | jerr << "Unable to get TOF_ADC_TSCALE from /TOF/digi_scales !" << endl; |
109 | } |
110 | |
111 | |
112 | map<string,double> base_time_offset; |
113 | if (eventLoop->GetCalib("/TOF/base_time_offset",base_time_offset)) |
114 | jout << "Error loading /TOF/base_time_offset !" << endl; |
115 | if (base_time_offset.find("TOF_BASE_TIME_OFFSET") != base_time_offset.end()) |
116 | t_base = base_time_offset["TOF_BASE_TIME_OFFSET"]; |
117 | else |
118 | jerr << "Unable to get TOF_BASE_TIME_OFFSET from /TOF/base_time_offset !" << endl; |
119 | |
120 | if (base_time_offset.find("TOF_TDC_BASE_TIME_OFFSET") != base_time_offset.end()) |
121 | t_base_tdc = base_time_offset["TOF_TDC_BASE_TIME_OFFSET"]; |
122 | else |
123 | jerr << "Unable to get TOF_TDC_BASE_TIME_OFFSET from /TOF/base_time_offset !" << endl; |
124 | |
125 | |
126 | if(eventLoop->GetCalib("TOF/pedestals", raw_adc_pedestals)) |
127 | jout << "Error loading /TOF/pedestals !" << endl; |
128 | if(eventLoop->GetCalib("TOF/gains", raw_adc_gains)) |
129 | jout << "Error loading /TOF/gains !" << endl; |
130 | if(eventLoop->GetCalib("TOF/adc_timing_offsets", raw_adc_offsets)) |
131 | jout << "Error loading /TOF/adc_timing_offsets !" << endl; |
132 | if(eventLoop->GetCalib("TOF/timing_offsets", raw_tdc_offsets)) |
133 | jout << "Error loading /TOF/timing_offsets !" << endl; |
134 | if(eventLoop->GetCalib("TOF/timewalk_parms", timewalk_parameters)) |
135 | jout << "Error loading /TOF/timewalk_parms !" << endl; |
136 | |
137 | FillCalibTable(adc_pedestals, raw_adc_pedestals, tofGeom); |
138 | FillCalibTable(adc_gains, raw_adc_gains, tofGeom); |
139 | FillCalibTable(adc_time_offsets, raw_adc_offsets, tofGeom); |
140 | FillCalibTable(tdc_time_offsets, raw_tdc_offsets, tofGeom); |
141 | |
142 | if(eventLoop->GetCalib("TOF/adc2E", raw_adc2E)) |
143 | jout << "Error loading /TOF/adc2E !" << endl; |
144 | |
145 | for (unsigned int n=0; n<raw_adc2E.size(); n++){ |
146 | adc2E[n] = raw_adc2E[n]; |
147 | } |
148 | |
149 | |
150 | |
151 | |
152 | |
153 | |
154 | |
155 | |
156 | return NOERROR; |
157 | } |
158 | |
159 | |
160 | |
161 | |
162 | jerror_t DTOFHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) |
163 | { |
164 | |
165 | |
166 | |
167 | |
168 | |
169 | |
170 | |
171 | |
172 | |
173 | const DTTabUtilities* locTTabUtilities = NULL__null; |
174 | loop->GetSingle(locTTabUtilities); |
175 | |
176 | |
177 | vector<const DTOFDigiHit*> digihits; |
178 | loop->Get(digihits); |
179 | for(unsigned int i=0; i<digihits.size(); i++){ |
180 | const DTOFDigiHit *digihit = digihits[i]; |
181 | |
182 | |
183 | if(digihit->datasource == 1) { |
184 | |
185 | |
186 | |
187 | const Df250PulsePedestal* PPobj = NULL__null; |
188 | digihit->GetSingle(PPobj); |
189 | if (PPobj != NULL__null) { |
190 | if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue; |
191 | } else |
192 | continue; |
193 | |
194 | |
195 | } |
196 | |
197 | if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF)) |
198 | continue; |
199 | |
200 | |
201 | |
202 | |
203 | double pedestal = GetConstant(adc_pedestals, digihit); |
| Value stored to 'pedestal' during its initialization is never read |
204 | double nsamples_integral = digihit->nsamples_integral; |
205 | double nsamples_pedestal = digihit->nsamples_pedestal; |
206 | |
207 | |
208 | if(nsamples_pedestal == 0) { |
209 | jerr << "DTOFDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl; |
210 | continue; |
211 | } |
212 | |
213 | if( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) { |
214 | pedestal = digihit->pedestal * (nsamples_integral/nsamples_pedestal); |
215 | } else { |
216 | continue; |
217 | } |
218 | |
219 | |
220 | |
221 | |
222 | double A = (double)digihit->pulse_integral; |
223 | double T = (double)digihit->pulse_time; |
224 | T = t_scale * T - GetConstant(adc_time_offsets, digihit) + t_base; |
225 | double dA = A - pedestal; |
226 | |
227 | if (dA<0) continue; |
228 | |
229 | DTOFHit *hit = new DTOFHit; |
230 | hit->plane = digihit->plane; |
231 | hit->bar = digihit->bar; |
232 | hit->end = digihit->end; |
233 | hit->dE=dA; |
234 | |
235 | if(COSMIC_DATA) |
236 | hit->dE = (A - 55*pedestal); |
237 | |
238 | hit->t_TDC=numeric_limits<double>::quiet_NaN(); |
239 | hit->t_fADC=T; |
240 | hit->t = hit->t_fADC; |
241 | |
242 | hit->has_fADC=true; |
243 | hit->has_TDC=false; |
244 | |
245 | |
246 | |
247 | |
248 | |
249 | |
250 | |
251 | |
252 | hit->AddAssociatedObject(digihit); |
253 | |
254 | _data.push_back(hit); |
255 | } |
256 | |
257 | |
258 | vector<const DTOFTDCDigiHit*> tdcdigihits; |
259 | loop->Get(tdcdigihits); |
260 | |
261 | |
262 | |
263 | |
264 | |
265 | for(unsigned int i=0; i<tdcdigihits.size(); i++) |
266 | { |
267 | const DTOFTDCDigiHit *digihit = tdcdigihits[i]; |
268 | |
269 | |
270 | double T = locTTabUtilities->Convert_DigiTimeToNs_CAEN1290TDC(digihit); |
271 | T += t_base_tdc - GetConstant(tdc_time_offsets, digihit) + tdc_adc_time_offset; |
272 | |
273 | |
274 | |
275 | |
276 | |
277 | |
278 | |
279 | |
280 | DTOFHit *hit = FindMatch(digihit->plane, digihit->bar, digihit->end, T); |
281 | |
282 | if(!hit){ |
283 | hit = new DTOFHit; |
284 | hit->plane = digihit->plane; |
285 | hit->bar = digihit->bar; |
286 | hit->end = digihit->end; |
287 | hit->dE = 0.0; |
288 | hit->t_fADC=numeric_limits<double>::quiet_NaN(); |
289 | hit->has_fADC=false; |
290 | |
291 | _data.push_back(hit); |
292 | } |
293 | hit->has_TDC=true; |
294 | hit->t_TDC=T; |
295 | |
296 | if (hit->dE>0.){ |
297 | |
298 | |
299 | int id=88*hit->plane+44*hit->end+hit->bar-1; |
300 | double A=hit->dE; |
301 | double C1=timewalk_parameters[id][1]; |
302 | double C2=timewalk_parameters[id][2]; |
303 | double A0=timewalk_parameters[id][3]; |
304 | T-=C1*(pow(A,C2)-pow(A0,C2)); |
305 | } |
306 | hit->t=T; |
307 | |
308 | hit->AddAssociatedObject(digihit); |
309 | } |
310 | |
311 | |
312 | for (unsigned int i=0;i<_data.size();i++){ |
313 | int id=88*_data[i]->plane + 44*_data[i]->end + _data[i]->bar-1; |
314 | _data[i]->dE *= adc2E[id]; |
315 | |
316 | } |
317 | |
318 | return NOERROR; |
319 | } |
320 | |
321 | |
322 | |
323 | |
324 | DTOFHit* DTOFHit_factory::FindMatch(int plane, int bar, int end, double T) |
325 | { |
326 | DTOFHit* best_match = NULL__null; |
327 | |
328 | |
329 | |
330 | for(unsigned int i=0; i<_data.size(); i++){ |
331 | DTOFHit *hit = _data[i]; |
332 | |
333 | if(!isfinite(hit->t_fADC)) continue; |
334 | if(hit->plane != plane) continue; |
335 | if(hit->bar != bar) continue; |
336 | if(hit->end != end) continue; |
337 | |
338 | |
339 | double delta_T = fabs(T - hit->t); |
340 | if(delta_T > DELTA_T_ADC_TDC_MAX) continue; |
341 | |
342 | return hit; |
343 | |
344 | |
345 | if(best_match != NULL__null) { |
346 | if(delta_T < fabs(best_match->t - T)) |
347 | best_match = hit; |
348 | } else { |
349 | best_match = hit; |
350 | } |
351 | |
352 | } |
353 | |
354 | return best_match; |
355 | } |
356 | |
357 | |
358 | |
359 | |
360 | jerror_t DTOFHit_factory::erun(void) |
361 | { |
362 | return NOERROR; |
363 | } |
364 | |
365 | |
366 | |
367 | |
368 | jerror_t DTOFHit_factory::fini(void) |
369 | { |
370 | return NOERROR; |
371 | } |
372 | |
373 | |
374 | |
375 | |
376 | |
377 | void DTOFHit_factory::FillCalibTable(tof_digi_constants_t &table, vector<double> &raw_table, |
378 | const DTOFGeometry &tofGeom) |
379 | { |
380 | char str[256]; |
381 | int channel = 0; |
382 | |
383 | |
384 | table.clear(); |
385 | |
386 | for(int plane=0; plane<tofGeom.NLAYERS; plane++) { |
387 | int plane_index=2*TOF_NUM_BARS*plane; |
388 | table.push_back( vector< pair<double,double> >(TOF_NUM_BARS) ); |
389 | for(int bar=0; bar<TOF_NUM_BARS; bar++) { |
390 | table[plane][bar] |
391 | = pair<double,double>(raw_table[plane_index+bar], |
392 | raw_table[plane_index+TOF_NUM_BARS+bar]); |
393 | channel+=2; |
394 | } |
395 | } |
396 | |
397 | |
398 | if(channel != TOF_MAX_CHANNELS) { |
399 | sprintf(str, "Not enough channels for TOF table! channel=%d (should be %d)", |
400 | channel, TOF_MAX_CHANNELS); |
401 | cerr << str << endl; |
402 | throw JException(str); |
403 | } |
404 | } |
405 | |
406 | |
407 | |
408 | |
409 | |
410 | |
411 | |
412 | |
413 | |
414 | |
415 | |
416 | |
417 | |
418 | const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table, |
419 | const int in_plane, const int in_bar, const int in_end) const |
420 | { |
421 | char str[256]; |
422 | |
423 | if( (in_plane < 0) || (in_plane > TOF_NUM_PLANES)) { |
424 | sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_plane, TOF_NUM_PLANES); |
425 | cerr << str << endl; |
426 | throw JException(str); |
427 | } |
428 | if( (in_bar <= 0) || (in_bar > TOF_NUM_BARS)) { |
429 | sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_bar, TOF_NUM_BARS); |
430 | cerr << str << endl; |
431 | throw JException(str); |
432 | } |
433 | if( (in_end != 0) && (in_end != 1) ) { |
434 | sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_end); |
435 | cerr << str << endl; |
436 | throw JException(str); |
437 | } |
438 | |
439 | |
440 | |
441 | if(in_end == 0) { |
442 | return the_table[in_plane][in_bar].first; |
443 | } else { |
444 | return the_table[in_plane][in_bar].second; |
445 | } |
446 | } |
447 | |
448 | const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table, |
449 | const DTOFHit *in_hit) const |
450 | { |
451 | char str[256]; |
452 | |
453 | if( (in_hit->plane < 0) || (in_hit->plane > TOF_NUM_PLANES)) { |
454 | sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_hit->plane, TOF_NUM_PLANES); |
455 | cerr << str << endl; |
456 | throw JException(str); |
457 | } |
458 | if( (in_hit->bar <= 0) || (in_hit->bar > TOF_NUM_BARS)) { |
459 | sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_hit->bar, TOF_NUM_BARS); |
460 | cerr << str << endl; |
461 | throw JException(str); |
462 | } |
463 | if( (in_hit->end != 0) && (in_hit->end != 1) ) { |
464 | sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_hit->end); |
465 | cerr << str << endl; |
466 | throw JException(str); |
467 | } |
468 | |
469 | |
470 | |
471 | if(in_hit->end == 0) { |
472 | return the_table[in_hit->plane][in_hit->bar-1].first; |
473 | } else { |
474 | return the_table[in_hit->plane][in_hit->bar-1].second; |
475 | } |
476 | } |
477 | |
478 | const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table, |
479 | const DTOFDigiHit *in_digihit) const |
480 | { |
481 | char str[256]; |
482 | |
483 | if( (in_digihit->plane < 0) || (in_digihit->plane > TOF_NUM_PLANES)) { |
484 | sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->plane, TOF_NUM_PLANES); |
485 | cerr << str << endl; |
486 | throw JException(str); |
487 | } |
488 | if( (in_digihit->bar <= 0) || (in_digihit->bar > TOF_NUM_BARS)) { |
489 | sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->bar, TOF_NUM_BARS); |
490 | cerr << str << endl; |
491 | throw JException(str); |
492 | } |
493 | if( (in_digihit->end != 0) && (in_digihit->end != 1) ) { |
494 | sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_digihit->end); |
495 | cerr << str << endl; |
496 | throw JException(str); |
497 | } |
498 | |
499 | |
500 | |
501 | if(in_digihit->end == 0) { |
502 | return the_table[in_digihit->plane][in_digihit->bar-1].first; |
503 | } else { |
504 | return the_table[in_digihit->plane][in_digihit->bar-1].second; |
505 | } |
506 | } |
507 | |
508 | const double DTOFHit_factory::GetConstant( const tof_digi_constants_t &the_table, |
509 | const DTOFTDCDigiHit *in_digihit) const |
510 | { |
511 | char str[256]; |
512 | |
513 | if( (in_digihit->plane < 0) || (in_digihit->plane > TOF_NUM_PLANES)) { |
514 | sprintf(str, "Bad module # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->plane, TOF_NUM_PLANES); |
515 | cerr << str << endl; |
516 | throw JException(str); |
517 | } |
518 | if( (in_digihit->bar <= 0) || (in_digihit->bar > TOF_NUM_BARS)) { |
519 | sprintf(str, "Bad layer # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->bar, TOF_NUM_BARS); |
520 | cerr << str << endl; |
521 | throw JException(str); |
522 | } |
523 | if( (in_digihit->end != 0) && (in_digihit->end != 1) ) { |
524 | sprintf(str, "Bad end # requested in DTOFHit_factory::GetConstant()! requested=%d , should be 0-1", in_digihit->end); |
525 | cerr << str << endl; |
526 | throw JException(str); |
527 | } |
528 | |
529 | |
530 | |
531 | if(in_digihit->end == 0) { |
532 | return the_table[in_digihit->plane][in_digihit->bar-1].first; |
533 | } else { |
534 | return the_table[in_digihit->plane][in_digihit->bar-1].second; |
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 | |
574 | |
575 | |
576 | |
577 | |
578 | |