1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | #include <iostream> |
13 | #include <iomanip> |
14 | #include <math.h> |
15 | using namespace std; |
16 | |
17 | #include "DPSPair_factory.h" |
18 | #include "DPSHit.h" |
19 | using namespace jana; |
20 | |
21 | inline bool SortByTimeDifference(const DPSPair* pair1, const DPSPair* pair2) |
22 | { |
23 | double tdiff1 = fabs(pair1->ee.first->t-pair1->ee.second->t); |
24 | double tdiff2 = fabs(pair2->ee.first->t-pair2->ee.second->t); |
25 | return (tdiff1<tdiff2); |
26 | } |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | jerror_t DPSPair_factory::init(void) |
34 | { |
35 | DELTA_T_CLUST_MAX = 10.0; |
36 | DELTA_T_PAIR_MAX = 10.0; |
37 | |
38 | gPARMS->SetDefaultParameter("PSPair:DELTA_T_CLUST_MAX",DELTA_T_CLUST_MAX, |
39 | "Maximum difference in ns between hits in a cluster" |
40 | " in left and right arm of fine PS"); |
41 | |
42 | gPARMS->SetDefaultParameter("PSPair:DELTA_T_PAIR_MAX",DELTA_T_PAIR_MAX, |
43 | "Maximum difference in ns between a pair of hits" |
44 | " in left and right arm of fine PS"); |
45 | |
46 | return NOERROR; |
47 | } |
48 | |
49 | |
50 | |
51 | |
52 | jerror_t DPSPair_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) |
53 | { |
54 | return NOERROR; |
55 | } |
56 | |
57 | |
58 | |
59 | |
60 | jerror_t DPSPair_factory::evnt(JEventLoop *loop, uint64_t eventnumber) |
61 | { |
62 | |
63 | |
64 | |
65 | vector<const DPSHit*> hits; |
66 | loop->Get(hits); |
67 | |
68 | int debug = 0; |
69 | |
70 | tiles_left.clear(); |
71 | tiles_right.clear(); |
72 | |
73 | clust_left.clear(); |
74 | clust_right.clear(); |
75 | |
76 | |
77 | for (unsigned int ii = 0; ii < hits.size(); ii++) { |
78 | if( hits[ii]->arm == 0 ) { |
79 | tile tmp; |
80 | tmp.column = hits[ii]->column; |
81 | tmp.energy = hits[ii]->E; |
82 | tmp.integral = hits[ii]->integral; |
83 | tmp.pulse_peak = hits[ii]->pulse_peak; |
84 | tmp.time = hits[ii]->t; |
85 | tmp.used = 0; |
86 | tiles_left.push_back(tmp); |
87 | } |
88 | if( hits[ii]->arm == 1 ) { |
89 | tile tmp; |
90 | tmp.column = hits[ii]->column; |
91 | tmp.energy = hits[ii]->E; |
92 | tmp.integral = hits[ii]->integral; |
93 | tmp.pulse_peak = hits[ii]->pulse_peak; |
94 | tmp.time = hits[ii]->t; |
95 | tmp.used = 0; |
96 | tiles_right.push_back(tmp); |
97 | } |
98 | } |
99 | |
100 | sort(tiles_left.begin(), tiles_left.end(),SortByTile); |
101 | sort(tiles_right.begin(),tiles_right.end(),SortByTile); |
102 | |
103 | int last_tile = -1; |
104 | |
105 | vector<int> my_index; |
106 | |
107 | |
108 | if(debug) |
109 | cout << " Tiles Left Side = " << tiles_left.size() << endl; |
110 | |
111 | if(tiles_left.size() > 0){ |
112 | |
113 | for (unsigned int ii = 0; ii < tiles_left.size(); ii++){ |
114 | |
115 | my_index.clear(); |
116 | |
117 | if( tiles_left[ii].used == 1) continue; |
118 | |
119 | tiles_left[ii].used = 1; |
120 | |
121 | last_tile = tiles_left[ii].column; |
122 | |
123 | my_index.push_back(ii); |
124 | |
125 | for (unsigned int jj = ii + 1; jj < tiles_left.size(); jj++){ |
126 | if(fabs(tiles_left[ii].time - tiles_left[jj].time) < DELTA_T_CLUST_MAX){ |
127 | if( std::abs(tiles_left[jj].column - last_tile) <= 1){ |
128 | |
129 | tiles_left[jj].used = 1; |
130 | last_tile = tiles_left[jj].column; |
131 | |
132 | my_index.push_back(jj); |
133 | } |
134 | } |
135 | } |
136 | |
137 | clust tmp; |
138 | tmp.energy = 1; |
139 | tmp.hit_index = my_index; |
140 | |
141 | clust_left.push_back(tmp); |
142 | } |
143 | } |
144 | |
145 | |
146 | |
147 | last_tile = -1; |
| Value stored to 'last_tile' is never read |
148 | |
149 | if(debug){ |
150 | cout << endl; |
151 | cout << " Tiles Right Side = " << tiles_right.size() << endl; |
152 | } |
153 | |
154 | if(tiles_right.size() > 0){ |
155 | |
156 | for (unsigned int ii = 0; ii < tiles_right.size(); ii++){ |
157 | |
158 | my_index.clear(); |
159 | |
160 | if( tiles_right[ii].used == 1) continue; |
161 | |
162 | tiles_right[ii].used = 1; |
163 | |
164 | last_tile = tiles_right[ii].column; |
165 | |
166 | my_index.push_back(ii); |
167 | |
168 | for (unsigned int jj = ii + 1; jj < tiles_right.size(); jj++){ |
169 | if(fabs(tiles_right[ii].time - tiles_right[jj].time) < DELTA_T_CLUST_MAX){ |
170 | if( std::abs(tiles_right[jj].column - last_tile) <= 1){ |
171 | |
172 | tiles_right[jj].used = 1; |
173 | last_tile = tiles_right[jj].column; |
174 | |
175 | my_index.push_back(jj); |
176 | } |
177 | } |
178 | } |
179 | |
180 | clust tmp; |
181 | tmp.energy = 1; |
182 | tmp.hit_index = my_index; |
183 | |
184 | clust_right.push_back(tmp); |
185 | } |
186 | } |
187 | |
188 | |
189 | |
190 | |
191 | double max_integral = 0.; |
192 | double max_peak = 0.; |
193 | double max_time = 0.; |
194 | int max_column = -1; |
195 | |
196 | if( clust_left.size() > 0){ |
197 | for(unsigned ii = 0; ii < clust_left.size(); ii++){ |
198 | unsigned int nhits = clust_left[ii].hit_index.size(); |
199 | if(nhits > 0){ |
200 | double norm = 0.; |
201 | double en_tmp = 0.; |
202 | double time_tmp = 0.; |
203 | |
204 | max_integral = 0.; |
205 | |
206 | for(unsigned jj = 0; jj < nhits; jj++){ |
207 | int index = clust_left[ii].hit_index[jj]; |
208 | en_tmp += tiles_left[index].energy*tiles_left[index].integral; |
209 | time_tmp += tiles_left[index].time*tiles_left[index].integral; |
210 | norm += tiles_left[index].integral; |
211 | if(tiles_left[index].integral > max_integral) { |
212 | max_integral = tiles_left[index].integral; |
213 | max_peak = tiles_left[index].pulse_peak; |
214 | max_time = tiles_left[index].time; |
215 | max_column = tiles_left[index].column; |
216 | } |
217 | } |
218 | if(norm > 0){ |
219 | clust_left[ii].column = max_column; |
220 | clust_left[ii].energy = en_tmp/norm; |
221 | clust_left[ii].time = time_tmp/norm; |
222 | clust_left[ii].integral = max_integral; |
223 | clust_left[ii].pulse_peak = max_peak; |
224 | clust_left[ii].time_tile = max_time; |
225 | } else { |
226 | clust_left[ii].column = 0; |
227 | clust_left[ii].energy = 0; |
228 | clust_left[ii].time = 0; |
229 | clust_left[ii].integral = 0; |
230 | clust_left[ii].pulse_peak = 0; |
231 | clust_left[ii].time_tile = 0; |
232 | } |
233 | } |
234 | } |
235 | } |
236 | |
237 | max_integral = 0.; |
238 | max_peak = 0.; |
239 | max_time = 0.; |
240 | max_column = -1; |
241 | |
242 | if( clust_right.size() > 0){ |
243 | for(unsigned ii = 0; ii < clust_right.size(); ii++){ |
244 | unsigned int nhits = clust_right[ii].hit_index.size(); |
245 | if(nhits > 0){ |
246 | double norm = 0.; |
247 | double en_tmp = 0.; |
248 | double time_tmp = 0.; |
249 | |
250 | |
251 | max_integral = 0.; |
252 | |
253 | for(unsigned jj = 0; jj < nhits; jj++){ |
254 | int index = clust_right[ii].hit_index[jj]; |
255 | en_tmp += tiles_right[index].energy*tiles_right[index].integral; |
256 | time_tmp += tiles_right[index].time*tiles_right[index].integral; |
257 | norm += tiles_right[index].integral; |
258 | if(tiles_right[index].integral > max_integral) { |
259 | max_integral = tiles_right[index].integral; |
260 | max_peak = tiles_right[index].pulse_peak; |
261 | max_time = tiles_right[index].time; |
262 | max_column = tiles_right[index].column; |
263 | } |
264 | } |
265 | if(norm > 0){ |
266 | clust_right[ii].column = max_column; |
267 | clust_right[ii].energy = en_tmp/norm; |
268 | clust_right[ii].time = time_tmp/norm; |
269 | clust_right[ii].integral = max_integral; |
270 | clust_right[ii].pulse_peak = max_peak; |
271 | clust_right[ii].time_tile = max_time; |
272 | } else { |
273 | clust_right[ii].column = 0; |
274 | clust_right[ii].energy = 0; |
275 | clust_right[ii].time = 0; |
276 | clust_right[ii].integral = 0; |
277 | clust_right[ii].pulse_peak = 0; |
278 | clust_right[ii].time_tile = 0; |
279 | } |
280 | } |
281 | } |
282 | } |
283 | |
284 | |
285 | if(debug) |
286 | cout << " Number of Left clusters found = " << clust_left.size() << endl; |
287 | |
288 | for(unsigned ii = 0; ii < clust_left.size(); ii++){ |
289 | |
290 | if(debug) |
291 | cout << " Number of tiles inside the cluster = " << |
292 | clust_left[ii].hit_index.size() << endl; |
293 | |
294 | for(unsigned jj = 0; jj < clust_left[ii].hit_index.size(); jj++){ |
295 | int index_l = clust_left[ii].hit_index[jj]; |
296 | if(debug) |
297 | cout << " Index = " << index_l << |
298 | " Tile = " << tiles_left[index_l].column << |
299 | " Energy = " << tiles_left[index_l].energy << |
300 | " Integral = " << tiles_left[index_l].integral << endl; |
301 | } |
302 | |
303 | if(debug) |
304 | cout << " Cluster ENERGY: " << ii << " " << clust_left[ii].energy << |
305 | " TIME = " << clust_left[ii].time << endl; |
306 | |
307 | } |
308 | |
309 | |
310 | if(debug){ |
311 | cout << endl; |
312 | cout << endl; |
313 | cout << " Number of Right clusters found = " << clust_right.size() << endl; |
314 | } |
315 | |
316 | |
317 | for(unsigned ii = 0; ii < clust_right.size(); ii++){ |
318 | |
319 | if(debug) |
320 | cout << " Number of tiles inside the cluster = " << |
321 | clust_right[ii].hit_index.size() << endl; |
322 | |
323 | for(unsigned jj = 0; jj < clust_right[ii].hit_index.size(); jj++){ |
324 | int index_r = clust_right[ii].hit_index[jj]; |
325 | if(debug) |
326 | cout << " Index = " << index_r << |
327 | " Tile = " << tiles_right[index_r].column << |
328 | " Energy = " << tiles_right[index_r].energy << |
329 | " Integral = " << tiles_right[index_r].integral << endl; |
330 | } |
331 | |
332 | if(debug) |
333 | cout << " Cluster ENERGY: " << ii << " " << clust_right[ii].energy << |
334 | " TIME = " << clust_right[ii].time << endl; |
335 | |
336 | } |
337 | |
338 | |
339 | if( (clust_left.size() > 0) && (clust_right.size() > 0)) { |
340 | for(unsigned ii = 0; ii < clust_left.size(); ii++){ |
341 | for(unsigned jj = 0; jj < clust_right.size(); jj++){ |
342 | |
343 | if(fabs(clust_left[ii].time - clust_right[jj].time) < DELTA_T_PAIR_MAX){ |
344 | |
345 | if(debug){ |
346 | cout << endl; |
347 | cout << " PAIR FOUND: " << " E = " << |
348 | clust_left[ii].energy + clust_right[jj].energy << endl; |
349 | cout << endl; |
350 | } |
351 | |
352 | |
353 | DPSPair *pair = new DPSPair; |
354 | |
355 | pair->SetPair(clust_left[ii].column, clust_left[ii].pulse_peak, clust_left[ii].integral, clust_left[ii].time_tile, |
356 | clust_left[ii].hit_index.size(), clust_left[ii].energy, clust_left[ii].time, |
357 | clust_right[jj].column, clust_right[jj].pulse_peak, clust_right[jj].integral, clust_right[jj].time_tile, |
358 | clust_right[jj].hit_index.size(), clust_right[jj].energy, clust_right[jj].time ); |
359 | |
360 | _data.push_back(pair); |
361 | |
362 | } |
363 | |
364 | |
365 | } |
366 | } |
367 | } |
368 | |
369 | sort(_data.begin(),_data.end(), SortByTimeDifference); |
370 | |
371 | return NOERROR; |
372 | } |
373 | |
374 | |
375 | |
376 | |
377 | jerror_t DPSPair_factory::erun(void) |
378 | { |
379 | return NOERROR; |
380 | } |
381 | |
382 | |
383 | |
384 | |
385 | jerror_t DPSPair_factory::fini(void) |
386 | { |
387 | return NOERROR; |
388 | } |
389 | |
390 | bool DPSPair_factory::SortByTile(const tile &tile1, const tile &tile2) |
391 | { |
392 | if(tile1.column == tile2.column) return tile1.time < tile2.time; |
393 | return (tile1.column < tile2.column); |
394 | } |