1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | #include <iostream> |
8 | #include <iomanip> |
9 | using namespace std; |
10 | |
11 | #include "DPSHit_factory.h" |
12 | #include <DAQ/Df250PulsePedestal.h> |
13 | #include <DAQ/Df250PulseIntegral.h> |
14 | #include <TTAB/DTTabUtilities.h> |
15 | using namespace jana; |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | jerror_t DPSHit_factory::init(void) |
22 | { |
23 | ADC_THRESHOLD = 500.0; |
24 | gPARMS->SetDefaultParameter("PSHit:ADC_THRESHOLD",ADC_THRESHOLD, |
25 | "pedestal-subtracted pulse integral threshold"); |
26 | |
27 | |
28 | a_scale = 0.0001; |
29 | t_scale = 0.0625; |
30 | t_base = 0.; |
31 | |
32 | CHECK_FADC_ERRORS = false; |
33 | gPARMS->SetDefaultParameter("PSHit:CHECK_FADC_ERRORS", CHECK_FADC_ERRORS, "Set to 1 to reject hits with fADC250 errors, ser to 0 to keep these hits"); |
34 | |
35 | return NOERROR; |
36 | } |
37 | |
38 | |
39 | |
40 | |
41 | jerror_t DPSHit_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 | if(print_messages) jout << "In DPSHit_factory, loading constants..." << endl; |
56 | |
57 | |
58 | vector<const DPSGeometry*> psGeomVect; |
59 | eventLoop->Get(psGeomVect); |
60 | if (psGeomVect.size() < 1) |
61 | return OBJECT_NOT_AVAILABLE; |
62 | const DPSGeometry& psGeom = *(psGeomVect[0]); |
63 | |
64 | |
65 | map<string,double> scale_factors; |
66 | if (eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/digi_scales", scale_factors)) |
67 | jout << "Error loading /PHOTON_BEAM/pair_spectrometer/digi_scales !" << endl; |
68 | if (scale_factors.find("PS_ADC_ASCALE") != scale_factors.end()) |
69 | a_scale = scale_factors["PS_ADC_ASCALE"]; |
70 | else |
71 | jerr << "Unable to get PS_ADC_ASCALE from /PHOTON_BEAM/pair_spectrometer/digi_scales !" |
72 | << endl; |
73 | if (scale_factors.find("PS_ADC_TSCALE") != scale_factors.end()) |
74 | t_scale = scale_factors["PS_ADC_TSCALE"]; |
75 | else |
76 | jerr << "Unable to get PS_ADC_TSCALE from /PHOTON_BEAM/pair_spectrometer/digi_scales !" |
77 | << endl; |
78 | |
79 | |
80 | map<string,double> base_time_offset; |
81 | if (eventLoop->GetCalib("/PHOTON_BEAM/pair_spectrometer/base_time_offset",base_time_offset)) |
82 | jout << "Error loading /PHOTON_BEAM/pair_spectrometer/base_time_offset !" << endl; |
83 | if (base_time_offset.find("PS_FINE_BASE_TIME_OFFSET") != base_time_offset.end()) |
84 | t_base = base_time_offset["PS_FINE_BASE_TIME_OFFSET"]; |
85 | else |
86 | jerr << "Unable to get PS_FINE_BASE_TIME_OFFSET from /PHOTON_BEAM/pair_spectrometer/base_time_offset !" << endl; |
87 | |
88 | FillCalibTable(adc_pedestals, "/PHOTON_BEAM/pair_spectrometer/fine/adc_pedestals", psGeom); |
89 | FillCalibTable(adc_gains, "/PHOTON_BEAM/pair_spectrometer/fine/adc_gain_factors", psGeom); |
90 | FillCalibTable(adc_time_offsets, "/PHOTON_BEAM/pair_spectrometer/fine/adc_timing_offsets", psGeom); |
91 | |
92 | return NOERROR; |
93 | } |
94 | |
95 | |
96 | |
97 | |
98 | jerror_t DPSHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) |
99 | { |
100 | |
101 | |
102 | |
103 | |
104 | |
105 | |
106 | |
107 | |
108 | |
109 | |
110 | vector<const DPSGeometry*> psGeomVect; |
111 | eventLoop->Get(psGeomVect); |
112 | if (psGeomVect.size() < 1) |
113 | return OBJECT_NOT_AVAILABLE; |
114 | const DPSGeometry& psGeom = *(psGeomVect[0]); |
115 | |
116 | const DTTabUtilities* locTTabUtilities = NULL__null; |
117 | loop->GetSingle(locTTabUtilities); |
118 | |
119 | |
120 | vector<const DPSDigiHit*> digihits; |
121 | loop->Get(digihits); |
122 | char str[256]; |
123 | |
124 | for (unsigned int i=0; i < digihits.size(); i++){ |
125 | const DPSDigiHit *digihit = digihits[i]; |
126 | |
127 | |
128 | if( (digihit->arm < 0) || (digihit->arm >= psGeom.NUM_ARMS) ) { |
129 | sprintf(str, "DPSDigiHit arm out of range! arm=%d (should be 0-%d)", |
130 | digihit->arm, psGeom.NUM_ARMS); |
131 | throw JException(str); |
132 | } |
133 | if( (digihit->column <= 0) || (digihit->column > psGeom.NUM_FINE_COLUMNS) ) { |
134 | sprintf(str, "DPSDigiHit column out of range! column=%d (should be 0-%d)", |
135 | digihit->column, psGeom.NUM_FINE_COLUMNS); |
136 | throw JException(str); |
137 | } |
138 | |
139 | |
140 | if(CHECK_FADC_ERRORS && !locTTabUtilities->CheckFADC250_NoErrors(digihit->QF)) |
141 | continue; |
142 | |
143 | if(digihit->datasource == 1) { |
144 | |
145 | |
146 | |
147 | |
148 | const Df250PulsePedestal* PPobj = NULL__null; |
149 | digihit->GetSingle(PPobj); |
150 | if (PPobj != NULL__null){ |
151 | if (PPobj->pedestal == 0 || PPobj->pulse_peak == 0) continue; |
152 | } |
153 | else continue; |
154 | } |
155 | |
156 | |
157 | |
158 | double pedestal = GetConstant(adc_pedestals,digihit,psGeom); |
| Value stored to 'pedestal' during its initialization is never read |
159 | double nsamples_integral = (double)digihit->nsamples_integral; |
160 | double nsamples_pedestal = (double)digihit->nsamples_pedestal; |
161 | |
162 | |
163 | if(nsamples_pedestal == 0) { |
164 | jerr << "DPSDigiHit with nsamples_pedestal == 0 ! Event = " << eventnumber << endl; |
165 | continue; |
166 | } |
167 | |
168 | if ( (digihit->pedestal>0) && locTTabUtilities->CheckFADC250_PedestalOK(digihit->QF) ) { |
169 | |
170 | |
171 | |
172 | double ped_sum = (double)digihit->pedestal; |
173 | pedestal = ped_sum * nsamples_integral/nsamples_pedestal; |
174 | } |
175 | else continue; |
176 | |
177 | double single_sample_ped = pedestal/nsamples_pedestal; |
178 | |
179 | |
180 | double P = digihit->pulse_peak - single_sample_ped; |
181 | double A = (double)digihit->pulse_integral; |
182 | A -= pedestal; |
183 | |
184 | if (A < ADC_THRESHOLD) continue; |
185 | double T = (double)digihit->pulse_time; |
186 | |
187 | DPSHit *hit = new DPSHit; |
188 | |
189 | |
190 | |
191 | hit->arm = digihit->arm; |
192 | hit->column = digihit->column; |
193 | hit->integral = A; |
194 | hit->pulse_peak = P; |
195 | hit->npix_fadc = A * a_scale * GetConstant(adc_gains, digihit, psGeom); |
196 | hit->t = t_scale * T - GetConstant(adc_time_offsets, digihit, psGeom) + t_base; |
197 | hit->E = 0.5*(psGeom.getElow(digihit->arm,digihit->column) + psGeom.getEhigh(digihit->arm,digihit->column)); |
198 | |
199 | hit->AddAssociatedObject(digihit); |
200 | |
201 | _data.push_back(hit); |
202 | } |
203 | |
204 | |
205 | return NOERROR; |
206 | } |
207 | |
208 | |
209 | |
210 | |
211 | jerror_t DPSHit_factory::erun(void) |
212 | { |
213 | return NOERROR; |
214 | } |
215 | |
216 | |
217 | |
218 | |
219 | jerror_t DPSHit_factory::fini(void) |
220 | { |
221 | return NOERROR; |
222 | } |
223 | |
224 | |
225 | |
226 | |
227 | void DPSHit_factory::FillCalibTable(ps_digi_constants_t &table, string table_name, |
228 | const DPSGeometry &psGeom) |
229 | { |
230 | char str[256]; |
231 | |
232 | |
233 | if(eventLoop->GetCalib(table_name, table)) |
234 | jout << "Error loading " + table_name + " !" << endl; |
235 | |
236 | |
237 | |
238 | if( (int)table.size() != psGeom.NUM_FINE_COLUMNS ) { |
239 | sprintf(str, "PS table loaded with wrong size! number of columns=%d (should be %d)", |
240 | (int)table.size(), psGeom.NUM_FINE_COLUMNS ); |
241 | cerr << str << endl; |
242 | throw JException(str); |
243 | } |
244 | |
245 | for( int column=0; column < psGeom.NUM_FINE_COLUMNS; column++) { |
246 | if( (int)table[column].size() != psGeom.NUM_ARMS ) { |
247 | sprintf(str, "PS table loaded with wrong size! column=%d number of arms=%d (should be %d)", |
248 | column, (int)table[column].size(), psGeom.NUM_ARMS ); |
249 | cerr << str << endl; |
250 | throw JException(str); |
251 | } |
252 | } |
253 | } |
254 | |
255 | |
256 | |
257 | |
258 | |
259 | |
260 | |
261 | |
262 | |
263 | |
264 | const double DPSHit_factory::GetConstant( const ps_digi_constants_t &the_table, |
265 | const DPSGeometry::Arm in_arm, const int in_column, |
266 | const DPSGeometry &psGeom ) const |
267 | { |
268 | char str[256]; |
269 | |
270 | if( (in_arm != DPSGeometry::kNorth) && (in_arm != DPSGeometry::kSouth)) { |
271 | sprintf(str, "Bad arm requested in DPSHit_factory::GetConstant()! requested=%d , should be 0-%d", |
272 | static_cast<int>(in_arm), static_cast<int>(DPSGeometry::kSouth)); |
273 | cerr << str << endl; |
274 | throw JException(str); |
275 | } |
276 | if( (in_column <= 0) || (in_column > psGeom.NUM_FINE_COLUMNS)) { |
277 | sprintf(str, "Bad column # requested in DPSHit_factory::GetConstant()! requested=%d , should be 1-%d", in_column, psGeom.NUM_FINE_COLUMNS); |
278 | cerr << str << endl; |
279 | throw JException(str); |
280 | } |
281 | |
282 | |
283 | |
284 | if(in_arm == DPSGeometry::kNorth) { |
285 | return the_table[in_column-1][in_arm]; |
286 | } else { |
287 | return the_table[in_column-1][in_arm]; |
288 | } |
289 | } |
290 | |
291 | const double DPSHit_factory::GetConstant( const ps_digi_constants_t &the_table, |
292 | const DPSHit *in_hit, const DPSGeometry &psGeom ) const |
293 | { |
294 | char str[256]; |
295 | |
296 | if( (in_hit->arm != DPSGeometry::kNorth) && (in_hit->arm != DPSGeometry::kSouth)) { |
297 | sprintf(str, "Bad arm requested in DPSHit_factory::GetConstant()! requested=%d , should be 0-%d", |
298 | static_cast<int>(in_hit->arm), static_cast<int>(DPSGeometry::kSouth)); |
299 | cerr << str << endl; |
300 | throw JException(str); |
301 | } |
302 | if( (in_hit->column <= 0) || (in_hit->column > psGeom.NUM_FINE_COLUMNS)) { |
303 | sprintf(str, "Bad column # requested in DPSHit_factory::GetConstant()! requested=%d , should be 1-%d", in_hit->column, psGeom.NUM_FINE_COLUMNS); |
304 | cerr << str << endl; |
305 | throw JException(str); |
306 | } |
307 | |
308 | |
309 | |
310 | if(in_hit->arm == DPSGeometry::kNorth) { |
311 | return the_table[in_hit->column-1][in_hit->arm]; |
312 | } else { |
313 | return the_table[in_hit->column-1][in_hit->arm]; |
314 | } |
315 | } |
316 | |
317 | const double DPSHit_factory::GetConstant( const ps_digi_constants_t &the_table, |
318 | const DPSDigiHit *in_digihit, const DPSGeometry &psGeom) const |
319 | { |
320 | char str[256]; |
321 | |
322 | if( (in_digihit->arm != DPSGeometry::kNorth) && (in_digihit->arm != DPSGeometry::kSouth)) { |
323 | sprintf(str, "Bad arm requested in DPSHit_factory::GetConstant()! requested=%d , should be 0-%d", |
324 | static_cast<int>(in_digihit->arm), static_cast<int>(DPSGeometry::kSouth)); |
325 | cerr << str << endl; |
326 | throw JException(str); |
327 | } |
328 | if( (in_digihit->column <= 0) || (in_digihit->column > psGeom.NUM_FINE_COLUMNS)) { |
329 | sprintf(str, "Bad column # requested in DPSHit_factory::GetConstant()! requested=%d , should be 1-%d", in_digihit->column, psGeom.NUM_FINE_COLUMNS); |
330 | cerr << str << endl; |
331 | throw JException(str); |
332 | } |
333 | |
334 | |
335 | |
336 | if(in_digihit->arm == DPSGeometry::kNorth) { |
337 | return the_table[in_digihit->column-1][in_digihit->arm]; |
338 | } else { |
339 | return the_table[in_digihit->column-1][in_digihit->arm]; |
340 | } |
341 | } |
342 | |