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