1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | #include <iostream> |
10 | #include <iomanip> |
11 | using namespace std; |
12 | |
13 | #include "DTAGMDigiHit.h" |
14 | #include "DTAGMTDCDigiHit.h" |
15 | #include "DTAGMGeometry.h" |
16 | #include "DTAGMHit_factory_Calib.h" |
17 | #include <DAQ/Df250PulseIntegral.h> |
18 | #include <DAQ/Df250PulsePedestal.h> |
19 | #include <DAQ/Df250Config.h> |
20 | |
21 | using namespace jana; |
22 | |
23 | const int DTAGMHit_factory_Calib::k_fiber_good; |
24 | const int DTAGMHit_factory_Calib::k_fiber_bad; |
25 | const int DTAGMHit_factory_Calib::k_fiber_noisy; |
26 | |
27 | |
28 | |
29 | |
30 | jerror_t DTAGMHit_factory_Calib::init(void) |
31 | { |
32 | DELTA_T_ADC_TDC_MAX = 10.0; |
33 | USE_ADC = 0; |
34 | CUT_FACTOR = 1; |
35 | gPARMS->SetDefaultParameter("TAGMHit:DELTA_T_ADC_TDC_MAX", DELTA_T_ADC_TDC_MAX, |
36 | "Maximum difference in ns between a (calibrated) fADC time and" |
37 | " F1TDC time for them to be matched in a single hit"); |
38 | gPARMS->SetDefaultParameter("TAGMHit:CUT_FACTOR", CUT_FACTOR, "TAGM pulse integral cut factor, 0 = no cut"); |
39 | gPARMS->SetDefaultParameter("TAGMHit:USE_ADC", USE_ADC, "Use ADC times in TAGM"); |
40 | |
41 | CHECK_FADC_ERRORS = false; |
42 | gPARMS->SetDefaultParameter("TAGMHit:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits"); |
43 | |
44 | |
45 | fadc_a_scale = 0; |
46 | fadc_t_scale = 0; |
47 | t_base = 0; |
48 | t_tdc_base=0; |
49 | |
50 | |
51 | for (int row = 0; row <= TAGM_MAX_ROW5; ++row) { |
52 | for (int col = 0; col <= TAGM_MAX_COLUMN102; ++col) { |
53 | fadc_gains[row][col] = 0; |
54 | fadc_pedestals[row][col] = 0; |
55 | fadc_time_offsets[row][col] = 0; |
56 | tdc_time_offsets[row][col] = 0; |
57 | fiber_quality[row][col] = 0; |
58 | } |
59 | } |
60 | |
61 | return NOERROR; |
62 | } |
63 | |
64 | |
65 | |
66 | |
67 | jerror_t DTAGMHit_factory_Calib::brun(jana::JEventLoop *eventLoop, int32_t runnumber) |
68 | { |
69 | |
70 | static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
71 | static set<int> runs_announced; |
72 | pthread_mutex_lock(&print_mutex); |
73 | bool print_messages = false; |
74 | if(runs_announced.find(runnumber) == runs_announced.end()){ |
75 | print_messages = true; |
76 | runs_announced.insert(runnumber); |
77 | } |
78 | pthread_mutex_unlock(&print_mutex); |
79 | |
80 | |
81 | fadc_a_scale = 1.1; |
82 | fadc_t_scale = 0.0625; |
83 | t_base = 0.; |
84 | |
85 | pthread_mutex_unlock(&print_mutex); |
86 | if(print_messages) jout << "In DTAGMHit_factory_Calib, loading constants..." << std::endl; |
87 | |
88 | |
89 | map<string,double> base_time_offset; |
90 | if (eventLoop->GetCalib("/PHOTON_BEAM/microscope/base_time_offset",base_time_offset)) |
91 | jout << "Error loading /PHOTON_BEAM/microscope/base_time_offset !" << endl; |
92 | if (base_time_offset.find("TAGM_BASE_TIME_OFFSET") != base_time_offset.end()) |
93 | t_base = base_time_offset["TAGM_BASE_TIME_OFFSET"]; |
94 | else |
95 | jerr << "Unable to get TAGM_BASE_TIME_OFFSET from /PHOTON_BEAM/microscope/base_time_offset !" << endl; |
96 | if (base_time_offset.find("TAGM_TDC_BASE_TIME_OFFSET") != base_time_offset.end()) |
97 | t_tdc_base = base_time_offset["TAGM_TDC_BASE_TIME_OFFSET"]; |
98 | else |
99 | jerr << "Unable to get TAGM_TDC_BASE_TIME_OFFSET from /PHOTON_BEAM/microscope/base_time_offset !" << endl; |
100 | |
101 | if (load_ccdb_constants("fadc_gains", "gain", fadc_gains) && |
102 | load_ccdb_constants("fadc_pedestals", "pedestal", fadc_pedestals) && |
103 | load_ccdb_constants("fadc_time_offsets", "offset", fadc_time_offsets) && |
104 | load_ccdb_constants("tdc_time_offsets", "offset", tdc_time_offsets) && |
105 | load_ccdb_constants("fiber_quality", "code", fiber_quality) && |
106 | load_ccdb_constants("tdc_timewalk_corrections", "c0", tw_c0) && |
107 | load_ccdb_constants("tdc_timewalk_corrections", "c1", tw_c1) && |
108 | load_ccdb_constants("tdc_timewalk_corrections", "c2", tw_c2) && |
109 | load_ccdb_constants("tdc_timewalk_corrections", "threshold", tw_c3) && |
110 | load_ccdb_constants("tdc_timewalk_corrections", "reference", ref) && |
111 | load_ccdb_constants("integral_cuts", "integral", int_cuts)) |
112 | { |
113 | return NOERROR; |
114 | } |
115 | return UNRECOVERABLE_ERROR; |
116 | } |
117 | |
118 | |
119 | |
120 | |
121 | jerror_t DTAGMHit_factory_Calib::evnt(JEventLoop *loop, uint64_t eventnumber) |
122 | { |
123 | |
124 | |
125 | |
126 | |
127 | |
128 | |
129 | |
130 | |
131 | |
132 | |
133 | vector<const DTAGMGeometry*> tagmGeomVect; |
134 | eventLoop->Get( tagmGeomVect ); |
135 | if (tagmGeomVect.size() < 1) |
136 | return OBJECT_NOT_AVAILABLE; |
137 | const DTAGMGeometry& tagmGeom = *(tagmGeomVect[0]); |
138 | |
139 | const DTTabUtilities* locTTabUtilities = nullptr; |
140 | loop->GetSingle(locTTabUtilities); |
141 | |
142 | |
143 | vector<const DTAGMDigiHit*> digihits; |
144 | loop->Get(digihits); |
145 | for (unsigned int i=0; i < digihits.size(); i++) { |
146 | const DTAGMDigiHit *digihit = digihits[i]; |
147 | |
148 | |
149 | if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF)) |
150 | continue; |
151 | |
152 | |
153 | |
154 | double pedestal = fadc_pedestals[digihit->row][digihit->column]; |
155 | double nsamples_integral = (double)digihit->nsamples_integral; |
156 | double nsamples_pedestal = (double)digihit->nsamples_pedestal; |
157 | |
158 | |
159 | if(nsamples_pedestal == 0) { |
160 | jerr << "DTAGMDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl; |
161 | continue; |
162 | } |
163 | |
164 | if ( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) { |
165 | |
166 | |
167 | |
168 | double ped_sum = (double)digihit->pedestal; |
169 | pedestal = ped_sum * nsamples_integral/nsamples_pedestal; |
170 | } |
171 | double single_sample_ped = pedestal/nsamples_pedestal; |
172 | |
173 | double pulse_peak = 0.0; |
174 | if(digihit->datasource == 1) { |
175 | |
176 | |
177 | |
178 | |
179 | const Df250PulsePedestal* PPobj = nullptr; |
180 | digihit->GetSingle(PPobj); |
181 | if (PPobj != nullptr){ |
182 | if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue; |
183 | pulse_peak = PPobj->pulse_peak - PPobj->pedestal; |
184 | } |
185 | |
186 | |
187 | |
188 | } else { |
189 | |
190 | pulse_peak = digihit->pulse_peak - single_sample_ped; |
| Value stored to 'pulse_peak' is never read |
191 | } |
192 | |
193 | |
194 | int quality = fiber_quality[digihit->row][digihit->column]; |
195 | if (quality == k_fiber_bad || quality == k_fiber_noisy) |
196 | continue; |
197 | |
198 | |
199 | double P = digihit->pulse_peak - pedestal/nsamples_pedestal; |
200 | |
201 | |
202 | double A = digihit->pulse_integral; |
203 | double T = digihit->pulse_time; |
204 | A -= pedestal; |
205 | int row = digihit->row; |
206 | int column = digihit->column; |
207 | if (A < CUT_FACTOR*int_cuts[row][column]) continue; |
208 | |
209 | DTAGMHit *hit = new DTAGMHit; |
210 | hit->row = row; |
211 | hit->column = column; |
212 | double Elow = tagmGeom.getElow(column); |
213 | double Ehigh = tagmGeom.getEhigh(column); |
214 | hit->E = (Elow + Ehigh)/2; |
215 | |
216 | hit->integral = A; |
217 | hit->pulse_peak = P; |
218 | hit->npix_fadc = A * fadc_a_scale * fadc_gains[row][column]; |
219 | hit->time_fadc = T * fadc_t_scale - fadc_time_offsets[row][column] + t_base; |
220 | hit->t=hit->time_fadc; |
221 | hit->has_TDC=false; |
222 | hit->has_fADC=true; |
223 | |
224 | hit->AddAssociatedObject(digihit); |
225 | _data.push_back(hit); |
226 | } |
227 | |
228 | |
229 | |
230 | |
231 | vector<const DTAGMTDCDigiHit*> tdcdigihits; |
232 | loop->Get(tdcdigihits); |
233 | for (unsigned int i=0; i < tdcdigihits.size(); i++) { |
234 | const DTAGMTDCDigiHit *digihit = tdcdigihits[i]; |
235 | |
236 | |
237 | int row = digihit->row; |
238 | int column = digihit->column; |
239 | double T = locTTabUtilities->Convert_DigiTimeToNs_F1TDC(digihit) - tdc_time_offsets[row][column] + t_tdc_base; |
240 | |
241 | |
242 | |
243 | DTAGMHit *hit = nullptr; |
244 | for (unsigned int j=0; j < _data.size(); ++j) { |
245 | if (_data[j]->row == row && _data[j]->column == column && |
246 | fabs(T - _data[j]->time_fadc) < DELTA_T_ADC_TDC_MAX) |
247 | { |
248 | hit = _data[j]; |
249 | } |
250 | } |
251 | if (hit == nullptr) { |
252 | hit = new DTAGMHit; |
253 | hit->row = row; |
254 | hit->column = column; |
255 | double Elow = tagmGeom.getElow(column); |
256 | double Ehigh = tagmGeom.getEhigh(column); |
257 | hit->E = (Elow + Ehigh)/2; |
258 | hit->time_fadc = 0; |
259 | hit->npix_fadc = 0; |
260 | hit->integral = 0; |
261 | hit->pulse_peak = 0; |
262 | hit->has_fADC=false; |
263 | _data.push_back(hit); |
264 | } |
265 | hit->time_tdc=T; |
266 | hit->has_TDC=true; |
267 | |
268 | |
269 | double P = hit->pulse_peak; |
270 | double c0 = tw_c0[row][column]; |
271 | double c1 = tw_c1[row][column]; |
272 | double c2 = tw_c2[row][column]; |
273 | |
274 | double c3 = tw_c3[row][column]; |
275 | double t0 = ref[row][column]; |
276 | |
277 | if (P > 0) { |
278 | |
279 | T -= c1*pow(1/(P+c3),c2) - (t0 - c0); |
280 | } |
281 | |
282 | |
283 | if (USE_ADC && hit->has_fADC) hit->t = hit->time_fadc; |
284 | else hit->t = T; |
285 | |
286 | hit->AddAssociatedObject(digihit); |
287 | } |
288 | |
289 | return NOERROR; |
290 | } |
291 | |
292 | |
293 | |
294 | |
295 | jerror_t DTAGMHit_factory_Calib::erun(void) |
296 | { |
297 | return NOERROR; |
298 | } |
299 | |
300 | |
301 | |
302 | |
303 | jerror_t DTAGMHit_factory_Calib::fini(void) |
304 | { |
305 | return NOERROR; |
306 | } |
307 | |
308 | |
309 | |
310 | |
311 | bool DTAGMHit_factory_Calib::load_ccdb_constants( |
312 | std::string table_name, |
313 | std::string column_name, |
314 | double result[TAGM_MAX_ROW5+1][TAGM_MAX_COLUMN102+1]) |
315 | { |
316 | std::vector< std::map<std::string, double> > table; |
317 | std::string ccdb_key = "/PHOTON_BEAM/microscope/" + table_name; |
318 | if (eventLoop->GetCalib(ccdb_key, table)) |
319 | { |
320 | jout << "Error loading " << ccdb_key << " from ccdb!" << std::endl; |
321 | return false; |
322 | } |
323 | for (unsigned int i=0; i < table.size(); ++i) { |
324 | int row = (table[i])["row"]; |
325 | int col = (table[i])["column"]; |
326 | result[row][col] = (table[i])[column_name]; |
327 | } |
328 | return true; |
329 | } |