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