1 | /** |
---|
2 | * |
---|
3 | * Author: John Idarraga <idarraga@cern.ch> , 2009 |
---|
4 | * Medipix Group, Universite de Montreal |
---|
5 | * |
---|
6 | */ |
---|
7 | |
---|
8 | #ifndef BlobsFinder_cpp |
---|
9 | #define BlobsFinder_cpp |
---|
10 | |
---|
11 | #include "BlobsFinder.h" |
---|
12 | #include "MAFTools/MAFTools.h" |
---|
13 | |
---|
14 | using namespace MSG; |
---|
15 | |
---|
16 | ClassImp(BlobsFinder) |
---|
17 | ClassImp(AllBlobsContainer) |
---|
18 | |
---|
19 | /* |
---|
20 | * AllBlobsContainer constructor |
---|
21 | * |
---|
22 | */ |
---|
23 | AllBlobsContainer::AllBlobsContainer(MediPixAlgo * algo) : |
---|
24 | CandidateContainer(algo) { |
---|
25 | |
---|
26 | } |
---|
27 | |
---|
28 | void AllBlobsContainer::SetBlobType(Int_t id, blobtype type){ |
---|
29 | allBlobs[id].SetBlobType(type); |
---|
30 | } |
---|
31 | |
---|
32 | blobtype AllBlobsContainer::GetBlobType(Int_t id) { |
---|
33 | return allBlobs[id].GetBlobType(); |
---|
34 | } |
---|
35 | |
---|
36 | TString AllBlobsContainer::GetTypeAsString(Int_t id) { |
---|
37 | return allBlobs[id].GetTypeAsString(); |
---|
38 | } |
---|
39 | |
---|
40 | void blob::SetBlobType(blobtype type) |
---|
41 | { |
---|
42 | btype = type; |
---|
43 | } |
---|
44 | |
---|
45 | /** |
---|
46 | * BlobsFinder constructor |
---|
47 | * |
---|
48 | */ |
---|
49 | BlobsFinder::BlobsFinder(){ |
---|
50 | |
---|
51 | aB = 0; |
---|
52 | // init parameters |
---|
53 | m_spiralEnds = NEIGHBORS_DISC_0_PIXELS; |
---|
54 | // generally safe to exclude the border line of pixels |
---|
55 | m_excludeBorder = 1; |
---|
56 | // minimum Chi2 required to perform rotation |
---|
57 | m_chisquare_OverDof_rotation = 1.0; |
---|
58 | // no discontinuity in blob's pixel members |
---|
59 | m_discontinuityTolerance = 0; |
---|
60 | |
---|
61 | // lvl1 if available |
---|
62 | m_lvl1Cut = -1; // -1: no lvl1 available |
---|
63 | // rejected by level1 |
---|
64 | m_bP_rejectedByLvl1 = 0; |
---|
65 | } |
---|
66 | |
---|
67 | /** |
---|
68 | * Calculate length of the spiral through the formula : |
---|
69 | * (2*tolerancePixels + 2)*4 + recursive_call(--tolerancePixels) |
---|
70 | * |
---|
71 | * disc --> length |
---|
72 | * 0 8 |
---|
73 | * 1 24 |
---|
74 | * 2 48 |
---|
75 | * 3 80 |
---|
76 | * and so on ... |
---|
77 | * |
---|
78 | */ |
---|
79 | Int_t BlobsFinder::GetLengthOfSpiral(Int_t tolerancePixels){ |
---|
80 | |
---|
81 | if(tolerancePixels == -1) |
---|
82 | return 0; |
---|
83 | else |
---|
84 | return (2*tolerancePixels + 2)*4 + |
---|
85 | GetLengthOfSpiral(--tolerancePixels); |
---|
86 | |
---|
87 | } |
---|
88 | |
---|
89 | /** |
---|
90 | * Set length of the spiral |
---|
91 | * |
---|
92 | */ |
---|
93 | void BlobsFinder::SetDiscontinuityTolerance(Int_t disc = 0){ |
---|
94 | |
---|
95 | m_discontinuityTolerance = disc; |
---|
96 | m_spiralEnds = GetLengthOfSpiral(disc); |
---|
97 | |
---|
98 | if(disc >= MAX_DISCONTINUITY) |
---|
99 | { |
---|
100 | Log << MSG::ERROR << "Trying to set a discontinuity value probaly too high ? --> " << disc << endreq; |
---|
101 | Log << MSG::ERROR << " check BlobsFinder::SetDiscontinuityTolerance." << endreq; |
---|
102 | exit(1); |
---|
103 | } |
---|
104 | |
---|
105 | } |
---|
106 | |
---|
107 | /** |
---|
108 | * Called at the beginning of the run, before the first frame is |
---|
109 | * processed. |
---|
110 | * |
---|
111 | */ |
---|
112 | void BlobsFinder::Init(){ |
---|
113 | |
---|
114 | // a few branches for this algorithm's Tree |
---|
115 | getMyTree()->Branch("nBlobs", &m_bP_nBlobs , "nBlobs/I"); |
---|
116 | getMyTree()->Branch("nRejectedLvl1", &m_bP_rejectedByLvl1, "nRejectedLvl1/I"); |
---|
117 | |
---|
118 | getMyTree()->Branch("nPixels", m_bP_nPixels, "nPixels[nBlobs]/I"); |
---|
119 | getMyTree()->Branch("totalCharge", m_bP_totalCharge, "totalCharge[nBlobs]/I"); |
---|
120 | getMyTree()->Branch("width_x", m_bP_width_x, "width_x[nBlobs]/I"); |
---|
121 | getMyTree()->Branch("width_y", m_bP_width_y, "width_y[nBlobs]/I"); |
---|
122 | getMyTree()->Branch("geoCenter_x", m_bP_geoCenter_x, "geoCenter_x[nBlobs]/F"); |
---|
123 | getMyTree()->Branch("geoCenter_y", m_bP_geoCenter_y, "geoCenter_y[nBlobs]/F"); |
---|
124 | getMyTree()->Branch("weightedCenter_x", m_bP_weightedCenter_x, "weightedCenter_x[nBlobs]/F"); |
---|
125 | getMyTree()->Branch("weightedCenter_y", m_bP_weightedCenter_y, "weightedCenter_y[nBlobs]/F"); |
---|
126 | getMyTree()->Branch("boxArea", m_bP_boxArea, "boxArea[nBlobs]/I"); |
---|
127 | getMyTree()->Branch("circleArea", m_bP_circleArea, "circleArea[nBlobs]/F"); |
---|
128 | getMyTree()->Branch("rotAngle", m_bP_rotAngle, "rotAngle[nBlobs]/F"); |
---|
129 | getMyTree()->Branch("ellipseArea", m_bP_ellipseArea, "ellipseArea[nBlobs]/F"); |
---|
130 | getMyTree()->Branch("balanceToMin", m_bP_balanceToMin, "balanceToMin[nBlobs]/F"); |
---|
131 | getMyTree()->Branch("balanceToMax", m_bP_balanceToMax, "balanceToMax[nBlobs]/F"); |
---|
132 | |
---|
133 | // register some values for manipulation from the MPXViewer |
---|
134 | RegisterConfigurationValue(&m_excludeBorder, "border"); |
---|
135 | RegisterConfigurationValue(&m_discontinuityTolerance, "discontinuity"); |
---|
136 | RegisterConfigurationValue(&m_lvl1Cut, "lvl1cut"); |
---|
137 | |
---|
138 | //RegisterConfigurationValue(&m_spiralEnds, "spiral"); |
---|
139 | |
---|
140 | } |
---|
141 | |
---|
142 | /** |
---|
143 | * Called once per frame |
---|
144 | * |
---|
145 | */ |
---|
146 | void BlobsFinder::Execute(){ |
---|
147 | |
---|
148 | // Get dimensions |
---|
149 | Int_t xDim = GetMatrixXdim(); |
---|
150 | Int_t yDim = GetMatrixYdim(); |
---|
151 | |
---|
152 | // Set limits |
---|
153 | limits m_padLimits; |
---|
154 | m_padLimits.colmin = m_excludeBorder; |
---|
155 | m_padLimits.colmax = xDim - m_excludeBorder; |
---|
156 | m_padLimits.rowmin = m_excludeBorder; |
---|
157 | m_padLimits.rowmax = yDim - m_excludeBorder; |
---|
158 | Int_t rowItr = 0; |
---|
159 | |
---|
160 | // Recalculate lenght or spiral in case discontinuity tolerance |
---|
161 | // changed. |
---|
162 | m_spiralEnds = GetLengthOfSpiral(m_discontinuityTolerance); |
---|
163 | |
---|
164 | // initialize One blob data |
---|
165 | ClearOneBlobData(); |
---|
166 | |
---|
167 | // new object containing all the blobs, ready to go to StoreGate when's ready |
---|
168 | aB = new AllBlobsContainer((MediPixAlgo *) this); |
---|
169 | Log << MSG::INFO << "container --> " << aB << endreq; |
---|
170 | |
---|
171 | // Loop over the whole matrix except the selected border |
---|
172 | for(Int_t colItr = m_padLimits.colmin ; colItr < m_padLimits.colmax ; colItr++) |
---|
173 | { |
---|
174 | for(rowItr = m_padLimits.rowmin ; rowItr < m_padLimits.rowmax ; rowItr++) |
---|
175 | { |
---|
176 | if(GetMatrixElement(colItr, rowItr) > 0) |
---|
177 | { |
---|
178 | if(StartBlob(colItr, rowItr)) |
---|
179 | { |
---|
180 | // calculate one blob's properties |
---|
181 | int retVal = CalculateProperties(); |
---|
182 | DumpProperties(); |
---|
183 | //DumpInfo(); |
---|
184 | // fill blobs vector ... this guy's going to StoreGate |
---|
185 | if(retVal == __BLOB_GOOD_TO_STORE) aB->allBlobs.push_back(oneBlob); |
---|
186 | |
---|
187 | // any rejection |
---|
188 | if(retVal == __REJECTED_BY_LEVEL_1) m_bP_rejectedByLvl1++; |
---|
189 | // clean up current blob |
---|
190 | ClearOneBlobData(); |
---|
191 | } |
---|
192 | } |
---|
193 | } |
---|
194 | } |
---|
195 | |
---|
196 | Log << MSG::INFO << aB->allBlobs.size() << " blobs found in frame " << GetFrameId() << endreq; |
---|
197 | |
---|
198 | // pull all blobs to store gage |
---|
199 | PullToStoreGateAccess(aB, MPXDefs::DO_NOT_SERIALIZE_ME); |
---|
200 | |
---|
201 | // Filling up information to put into this algo's ntuple. |
---|
202 | m_bP_nBlobs = aB->allBlobs.size(); |
---|
203 | if(getMyTree()->GetNbranches() > 0) CopyPropertiesToNtupleData(); |
---|
204 | |
---|
205 | // fill ntuple |
---|
206 | getMyTree()->Fill(); |
---|
207 | |
---|
208 | // clean up all blobs |
---|
209 | //aB->allBlobs.clear(); |
---|
210 | |
---|
211 | // rewind |
---|
212 | m_bP_nBlobs = 0; |
---|
213 | |
---|
214 | // delete the blobs container |
---|
215 | //delete aB; |
---|
216 | |
---|
217 | } |
---|
218 | |
---|
219 | /** |
---|
220 | * Called at the end of the run, after the last frame is processed. |
---|
221 | * |
---|
222 | */ |
---|
223 | void BlobsFinder::Finalize(){ |
---|
224 | |
---|
225 | } |
---|
226 | |
---|
227 | /** |
---|
228 | * Single blobs are treated as a transient structure. Here I clean |
---|
229 | * up that information after the blob was copied somewhere else. |
---|
230 | * |
---|
231 | */ |
---|
232 | void BlobsFinder::ClearOneBlobData(){ |
---|
233 | |
---|
234 | // clean up |
---|
235 | oneBlob.blobContent.clear(); |
---|
236 | oneBlob.blobContentSet.clear(); |
---|
237 | |
---|
238 | oneBlob.bP.nPixels = 0; |
---|
239 | oneBlob.bP.totalCharge = 0; |
---|
240 | oneBlob.bP.width_x = 0; |
---|
241 | oneBlob.bP.width_y = 0; |
---|
242 | oneBlob.bP.geoCenter_x = 0; |
---|
243 | oneBlob.bP.geoCenter_y = 0; |
---|
244 | oneBlob.bP.weightedCenter_x = 0; |
---|
245 | oneBlob.bP.weightedCenter_y = 0; |
---|
246 | oneBlob.bP.boxArea = 0; |
---|
247 | oneBlob.bP.circleArea = 0.; |
---|
248 | oneBlob.bP.minToGeoCenter = 0.; |
---|
249 | oneBlob.bP.maxToGeoCenter = 0.; |
---|
250 | |
---|
251 | oneBlob.bP.lvl1.clear(); |
---|
252 | |
---|
253 | } |
---|
254 | |
---|
255 | /** |
---|
256 | * Start point of a blob |
---|
257 | * |
---|
258 | */ |
---|
259 | bool BlobsFinder::StartBlob(Int_t x_ini, Int_t y_ini){ |
---|
260 | |
---|
261 | // start point of a blob |
---|
262 | pair<Int_t, Int_t> startPoint = make_pair(x_ini, y_ini); |
---|
263 | |
---|
264 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
265 | // Check if a blob containing this point has already been saved |
---|
266 | // I only need to find the 'startPoint' in 'bloblContentSet' in each blob |
---|
267 | vector<blob>::iterator blobsItr = aB->allBlobs.begin(); |
---|
268 | set< pair<Int_t, Int_t> >::iterator itr; |
---|
269 | for( ; blobsItr != aB->allBlobs.end() ; blobsItr++) |
---|
270 | { |
---|
271 | itr = (*blobsItr).blobContentSet.find(startPoint); |
---|
272 | if(itr != (*blobsItr).blobContentSet.end()) // point found in previous blob !!! |
---|
273 | return false; |
---|
274 | } |
---|
275 | |
---|
276 | // If we hit this point we are ready to start a new blob |
---|
277 | // Insert in <set>, push_back in <list>, first time |
---|
278 | oneBlob.blobContentSet.insert( startPoint ); |
---|
279 | oneBlob.blobContent.push_back( pair< pair<Int_t, Int_t>, Int_t > (startPoint, __CENTER)); |
---|
280 | oneBlob.blobItr = oneBlob.blobContent.begin(); |
---|
281 | |
---|
282 | // Follow |
---|
283 | FollowBlob(); |
---|
284 | |
---|
285 | return true; |
---|
286 | } |
---|
287 | |
---|
288 | /** |
---|
289 | * Using recursion to follow a blob up to its end. |
---|
290 | * |
---|
291 | */ |
---|
292 | bool BlobsFinder::FollowBlob(){ |
---|
293 | |
---|
294 | // test insertion in the set |
---|
295 | pair< set < pair<Int_t, Int_t> >::iterator , bool> rtest; |
---|
296 | |
---|
297 | // spiral starts, i.e. search through the neighbors |
---|
298 | RewindSpiral(); |
---|
299 | |
---|
300 | // spiral starts at |
---|
301 | pair<Int_t, Int_t> nextNeighbor = (*oneBlob.blobItr).first; |
---|
302 | |
---|
303 | while((*oneBlob.blobItr).second < m_spiralEnds) { |
---|
304 | |
---|
305 | // try current position, next neighbor |
---|
306 | //nextNeighbor = GetNextPosition( (*oneBlob.blobItr).first ); |
---|
307 | nextNeighbor = GetNextPosition( nextNeighbor ); |
---|
308 | |
---|
309 | // see if the next position does not fall in the pad |
---|
310 | if(IsUnsafePosition(nextNeighbor)) { |
---|
311 | // increment neighbor counter, and skip the rest withing the while loop |
---|
312 | (*oneBlob.blobItr).second++; |
---|
313 | continue; |
---|
314 | } |
---|
315 | |
---|
316 | |
---|
317 | //Log << MSG::LOOP_DEBUG << "next: " << nextNeighbor.first << " , " << nextNeighbor.second << " ---> " |
---|
318 | //<< GetMatrixElement(nextNeighbor) << endreq; |
---|
319 | |
---|
320 | // if there is something |
---|
321 | if(GetMatrixElement(nextNeighbor) > 0) |
---|
322 | { |
---|
323 | Log << MSG::LOOP_DEBUG << "next: " << nextNeighbor.first << " , " << nextNeighbor.second << " ---> " |
---|
324 | << GetMatrixElement(nextNeighbor) << endreq; |
---|
325 | |
---|
326 | // check the Set see if it existed |
---|
327 | rtest = oneBlob.blobContentSet.insert( nextNeighbor ); |
---|
328 | // if insertion is possible it means this point didn't exist. So, push_back new entry in the list |
---|
329 | if(rtest.second == true) { |
---|
330 | oneBlob.blobContent.push_back( pair< pair<Int_t, Int_t>, Int_t > (nextNeighbor, __CENTER) ); |
---|
331 | } |
---|
332 | |
---|
333 | } |
---|
334 | |
---|
335 | // increment neighbor counter |
---|
336 | (*oneBlob.blobItr).second++; |
---|
337 | } |
---|
338 | |
---|
339 | // go to the next in the list, if there are any left |
---|
340 | if(++oneBlob.blobItr != oneBlob.blobContent.end()) { |
---|
341 | FollowBlob(); |
---|
342 | } |
---|
343 | |
---|
344 | return true; |
---|
345 | } |
---|
346 | |
---|
347 | /** |
---|
348 | * Check if the position we are about to step in is safe, i.e. see if |
---|
349 | * it falls inside the pad. |
---|
350 | * |
---|
351 | */ |
---|
352 | bool BlobsFinder::IsUnsafePosition(pair<Int_t, Int_t> pos){ |
---|
353 | |
---|
354 | if(pos.first < 0 || pos.second < 0) |
---|
355 | return true; |
---|
356 | if(pos.first > GetMatrixXdim()-1 || pos.second > GetMatrixYdim()-1) |
---|
357 | return true; |
---|
358 | |
---|
359 | return false; |
---|
360 | } |
---|
361 | |
---|
362 | /** |
---|
363 | * Simple blob's contents dump |
---|
364 | * |
---|
365 | */ |
---|
366 | void BlobsFinder::DumpInfo(){ |
---|
367 | |
---|
368 | list< pair < pair<Int_t, Int_t>, Int_t > >::iterator blobItrLocal = oneBlob.blobContent.begin(); |
---|
369 | Int_t totalPixels = 0; |
---|
370 | for( ; blobItrLocal != oneBlob.blobContent.end() ; blobItrLocal++) |
---|
371 | { |
---|
372 | Log << MSG::DEBUG << (*blobItrLocal).first.first << " , " << (*blobItrLocal).first.second |
---|
373 | << " -- spiral scan --> " << (*blobItrLocal).second << " " << endreq; |
---|
374 | totalPixels++; |
---|
375 | } |
---|
376 | |
---|
377 | Log << MSG::DEBUG << "---------- " << totalPixels << " pixels in this blob" << endreq; |
---|
378 | |
---|
379 | } |
---|
380 | |
---|
381 | void BlobsFinder::DumpProperties(){ |
---|
382 | |
---|
383 | Log << MSG::DEBUG << "center : " << oneBlob.bP.geoCenter_x << " , " |
---|
384 | << oneBlob.bP.geoCenter_y << endreq; |
---|
385 | Log << MSG::DEBUG << "boxArea : " << oneBlob.bP.boxArea << endreq; |
---|
386 | Log << MSG::DEBUG << "circle : " << oneBlob.bP.circleArea << endreq; |
---|
387 | Log << MSG::DEBUG << "balance max : " << oneBlob.bP.balanceToMax << endreq; |
---|
388 | Log << MSG::DEBUG << "balance min : " << oneBlob.bP.balanceToMin << endreq; |
---|
389 | Log << MSG::DEBUG << "---------- " << oneBlob.bP.nPixels << " pixels in this blob" << endreq; |
---|
390 | |
---|
391 | } |
---|
392 | |
---|
393 | /** |
---|
394 | * When done looking around a pixel. Rewind spiral steering |
---|
395 | * information. |
---|
396 | * |
---|
397 | */ |
---|
398 | void BlobsFinder::RewindSpiral(){ |
---|
399 | |
---|
400 | m_spiral.xySwitch = X_STEP; |
---|
401 | m_spiral.localMax = 1; |
---|
402 | m_spiral.xCntr = 1; |
---|
403 | m_spiral.yCntr = 1; |
---|
404 | m_spiral.dir = 1; |
---|
405 | |
---|
406 | } |
---|
407 | |
---|
408 | /** |
---|
409 | * Follows a spiral shape around a pixel |
---|
410 | * |
---|
411 | */ |
---|
412 | pair<Int_t, Int_t> BlobsFinder::GetNextPosition(pair<Int_t, Int_t> current) { |
---|
413 | |
---|
414 | pair<Int_t, Int_t> newPos; |
---|
415 | |
---|
416 | if(m_spiral.xySwitch == X_STEP) |
---|
417 | { |
---|
418 | newPos.first = current.first - m_spiral.dir; |
---|
419 | newPos.second = current.second; |
---|
420 | m_spiral.xCntr++; |
---|
421 | |
---|
422 | if(m_spiral.xCntr > m_spiral.localMax) |
---|
423 | m_spiral.xySwitch = Y_STEP; |
---|
424 | } |
---|
425 | else if(m_spiral.xySwitch == Y_STEP) |
---|
426 | { |
---|
427 | newPos.second = current.second - m_spiral.dir; |
---|
428 | newPos.first = current.first; |
---|
429 | m_spiral.yCntr++; |
---|
430 | |
---|
431 | if(m_spiral.yCntr > m_spiral.localMax) |
---|
432 | { |
---|
433 | m_spiral.xySwitch = X_STEP; |
---|
434 | m_spiral.localMax++; |
---|
435 | m_spiral.xCntr = 1; |
---|
436 | m_spiral.yCntr = 1; |
---|
437 | m_spiral.dir *= -1; // change direction |
---|
438 | } |
---|
439 | } |
---|
440 | |
---|
441 | return newPos; |
---|
442 | } |
---|
443 | |
---|
444 | /** |
---|
445 | * Set variable m_excludeBorder. m_excludeBorder is the number of |
---|
446 | * pixels that I will ignore in the border of the pixel pad. |
---|
447 | * |
---|
448 | */ |
---|
449 | void BlobsFinder::SetBorderExclusion(Int_t m_ex_i = 1){ |
---|
450 | |
---|
451 | /* |
---|
452 | if(m_ex_i > GetMatrixXdim()/2 - 1) { |
---|
453 | Log << MSG::ERROR << "Trying to set a border too big ... abort" << endreq; |
---|
454 | exit(1); |
---|
455 | } |
---|
456 | */ |
---|
457 | m_excludeBorder = m_ex_i; |
---|
458 | |
---|
459 | } |
---|
460 | |
---|
461 | /** |
---|
462 | * Calculate properties of one blob once it is complete |
---|
463 | * |
---|
464 | */ |
---|
465 | int BlobsFinder::CalculateProperties(){ |
---|
466 | |
---|
467 | list< pair < pair<Int_t, Int_t>, Int_t > >::iterator listItr; |
---|
468 | |
---|
469 | // max an min, x and y in the pad |
---|
470 | Int_t max_x = 0; |
---|
471 | Int_t min_x = GetMatrixXdim()-1; |
---|
472 | Int_t max_y = 0; |
---|
473 | Int_t min_y = GetMatrixYdim()-1; |
---|
474 | |
---|
475 | // number of pixels in this blob |
---|
476 | oneBlob.bP.nPixels = oneBlob.blobContent.size(); |
---|
477 | |
---|
478 | // calculate total charge |
---|
479 | oneBlob.bP.totalCharge = 0; |
---|
480 | |
---|
481 | Int_t xcoor = -1; |
---|
482 | Int_t ycoor = -1; |
---|
483 | set< pair<Int_t, Int_t> > blobSet = oneBlob.GetContentSet(); |
---|
484 | oneBlob.bP.nInnerPixels = 0; |
---|
485 | Int_t inthits = 0; |
---|
486 | |
---|
487 | for(listItr = oneBlob.blobContent.begin() ; listItr != oneBlob.blobContent.end() ; listItr++) |
---|
488 | { |
---|
489 | // <x, y> |
---|
490 | oneBlob.bP.totalCharge += GetMatrixElement((*listItr).first); |
---|
491 | xcoor = (*listItr).first.first; |
---|
492 | ycoor = (*listItr).first.second; |
---|
493 | |
---|
494 | if(xcoor < min_x) min_x = xcoor; |
---|
495 | if(xcoor > max_x) max_x = xcoor; |
---|
496 | if(ycoor < min_y) min_y = ycoor; |
---|
497 | if(ycoor > max_y) max_y = ycoor; |
---|
498 | |
---|
499 | // Get here the number of inner pixels, cross sape. |
---|
500 | // A pixel is said to be an inner pixel if it is found |
---|
501 | // in the following configuration. |
---|
502 | // * |
---|
503 | // * * * |
---|
504 | // * |
---|
505 | set< pair<Int_t, Int_t> >::iterator itri = blobSet.begin(); |
---|
506 | inthits = 0; |
---|
507 | for( ; itri != blobSet.end() ; itri++){ |
---|
508 | if((*itri).first == xcoor && (*itri).second == ycoor+1) inthits++; |
---|
509 | if((*itri).first == xcoor && (*itri).second == ycoor-1) inthits++; |
---|
510 | if((*itri).first == xcoor-1 && (*itri).second == ycoor ) inthits++; |
---|
511 | if((*itri).first == xcoor+1 && (*itri).second == ycoor ) inthits++; |
---|
512 | } |
---|
513 | |
---|
514 | if(inthits >= MIN_INNERPIXEL){ // internal pixel |
---|
515 | oneBlob.bP.nInnerPixels++; |
---|
516 | } |
---|
517 | } |
---|
518 | |
---|
519 | // widths of the blob |
---|
520 | oneBlob.bP.width_x = (max_x - min_x) + 1; // it is 1 if the width is 1 pixel |
---|
521 | oneBlob.bP.width_y = (max_y - min_y) + 1; // |
---|
522 | |
---|
523 | // geometrical center |
---|
524 | oneBlob.bP.geoCenter_x = (max_x + min_x + 1)/2.0;//TMath::CeilNint((max_x + min_x)/2.0); |
---|
525 | oneBlob.bP.geoCenter_y = (max_y + min_y + 1)/2.0;//TMath::CeilNint((max_y + min_y)/2.0); |
---|
526 | |
---|
527 | // calc min and max distance between all pixels and geoCenter |
---|
528 | CalcMinMaxGeoCenter(oneBlob.bP.geoCenter_x, oneBlob.bP.geoCenter_y); |
---|
529 | |
---|
530 | // Weighted center |
---|
531 | // works like the mean of an histogram |
---|
532 | // N = Sum_{i=0}^{M-1} H_{i} // where M is the total number of entries |
---|
533 | // // H_{i} is the histogram entry |
---|
534 | // The Mean: mu = 1/N Sum_{i=0}^{M-1} i H_{i} |
---|
535 | Float_t Nw = 0.; |
---|
536 | Float_t elem = 0.; |
---|
537 | for(listItr = oneBlob.blobContent.begin() ; listItr != oneBlob.blobContent.end() ; listItr++) |
---|
538 | { |
---|
539 | // <x, y> |
---|
540 | xcoor = (*listItr).first.first; |
---|
541 | ycoor = (*listItr).first.second; |
---|
542 | // calc N |
---|
543 | elem = (Float_t)GetMatrixElement((*listItr).first); |
---|
544 | Nw += elem; |
---|
545 | // calc mu |
---|
546 | oneBlob.bP.weightedCenter_x += (Float_t)xcoor * elem; |
---|
547 | oneBlob.bP.weightedCenter_y += (Float_t)ycoor * elem; |
---|
548 | |
---|
549 | // the lvl1 |
---|
550 | Int_t lvl1 = GetLVL1((*listItr).first); |
---|
551 | if(m_lvl1Cut != -1 && lvl1 >= m_lvl1Cut) { // reject this blob |
---|
552 | return __REJECTED_BY_LEVEL_1; |
---|
553 | } |
---|
554 | oneBlob.bP.lvl1.push_back(GetLVL1((*listItr).first)); |
---|
555 | |
---|
556 | } |
---|
557 | oneBlob.bP.weightedCenter_x /= Nw; |
---|
558 | oneBlob.bP.weightedCenter_y /= Nw; |
---|
559 | |
---|
560 | // calc box area |
---|
561 | Int_t wx = oneBlob.bP.width_x; |
---|
562 | Int_t wy = oneBlob.bP.width_y; |
---|
563 | oneBlob.bP.boxArea = wx*wy; |
---|
564 | |
---|
565 | |
---|
566 | /* |
---|
567 | * Some characteristics are calculated on a rotated blob |
---|
568 | * calculated only if the blob is not a Single, double, triple or quad hit. |
---|
569 | * Also some characteristics don't apply to these types. |
---|
570 | */ |
---|
571 | if(SingleHit() || DoubleHit() || TripleHit() || QuadHit()) |
---|
572 | { |
---|
573 | oneBlob.bP.circleArea = -1.; |
---|
574 | oneBlob.bP.ellipseArea = -1.; |
---|
575 | oneBlob.bP.ellipseA = -1.; |
---|
576 | oneBlob.bP.ellipseB = -1.; |
---|
577 | oneBlob.bP.rotAngle = -1.; |
---|
578 | oneBlob.bP.chisquare_OverDof = -1; |
---|
579 | oneBlob.bP.fitSlope = -1; |
---|
580 | oneBlob.bP.fitCut = -1; |
---|
581 | oneBlob.bP.balanceToMin = -1; |
---|
582 | oneBlob.bP.balanceToMax = -1; |
---|
583 | return __BLOB_GOOD_TO_STORE; |
---|
584 | } |
---|
585 | |
---|
586 | // CalculateBalance ... in this algorithm I duplicate the grid to avoid |
---|
587 | // meaningless information from certain structures. |
---|
588 | // I also get a guess for the linear regression start parameters. |
---|
589 | Float_t balanceToMin = 0., balanceToMax = 0., |
---|
590 | slopeGuess = 0., intersectionyGuess = 0., |
---|
591 | chisquare_OverDof = 0.; |
---|
592 | Int_t minQId = 0, maxQId = 0; |
---|
593 | GetBalance(min_x, max_x, min_y, max_y, |
---|
594 | & balanceToMin, & balanceToMax, |
---|
595 | & minQId, & maxQId, |
---|
596 | & slopeGuess, & intersectionyGuess); |
---|
597 | |
---|
598 | oneBlob.bP.balanceToMin = balanceToMin; |
---|
599 | oneBlob.bP.balanceToMax = balanceToMax; |
---|
600 | |
---|
601 | // Calculate linear regression |
---|
602 | Float_t fitSlope = 0.; |
---|
603 | Float_t fitCut = 0.; |
---|
604 | |
---|
605 | MAFTools::LinearRegression(&fitSlope, &fitCut, &chisquare_OverDof |
---|
606 | ,slopeGuess, intersectionyGuess |
---|
607 | ,oneBlob.GetContentSet()); |
---|
608 | |
---|
609 | oneBlob.bP.chisquare_OverDof = chisquare_OverDof; |
---|
610 | oneBlob.bP.rotAngle = -1.*TMath::ATan(fitSlope); |
---|
611 | oneBlob.bP.fitSlope = fitSlope; |
---|
612 | oneBlob.bP.fitCut = fitCut; |
---|
613 | |
---|
614 | // check if the object needs rotation |
---|
615 | if(oneBlob.bP.chisquare_OverDof > m_chisquare_OverDof_rotation) |
---|
616 | { |
---|
617 | |
---|
618 | // now rotate every point in the blob and find new max and min |
---|
619 | Int_t rot_max_x = -2*GetMatrixXdim()-1; |
---|
620 | Int_t rot_min_x = 2*GetMatrixXdim()-1; |
---|
621 | Int_t rot_max_y = -2*GetMatrixYdim()-1; |
---|
622 | Int_t rot_min_y = 2*GetMatrixYdim()-1; |
---|
623 | |
---|
624 | pair<Int_t, Int_t> rotPair; |
---|
625 | |
---|
626 | //Int_t totalItr = 0; |
---|
627 | for(listItr = oneBlob.blobContent.begin() ; listItr != oneBlob.blobContent.end() ; listItr++) |
---|
628 | { |
---|
629 | // <rot_x, rot_y> |
---|
630 | rotPair = RotateVector2D((*listItr).first.first, (*listItr).first.second, oneBlob.bP.rotAngle); |
---|
631 | if(rotPair.first < rot_min_x) rot_min_x = rotPair.first; |
---|
632 | if(rotPair.first > rot_max_x) rot_max_x = rotPair.first; |
---|
633 | if(rotPair.second < rot_min_y) rot_min_y = rotPair.second; |
---|
634 | if(rotPair.second > rot_max_y) rot_max_y = rotPair.second; |
---|
635 | |
---|
636 | //totalItr++; |
---|
637 | } |
---|
638 | |
---|
639 | Int_t rot_wx = (rot_max_x - rot_min_x) + 1; |
---|
640 | Int_t rot_wy = (rot_max_y - rot_min_y) + 1; |
---|
641 | |
---|
642 | oneBlob.bP.ellipseA = rot_wx; |
---|
643 | oneBlob.bP.ellipseB = rot_wy; |
---|
644 | |
---|
645 | // wx/2 is a, wy/2 is b. AreaElipse = Pi*a*b |
---|
646 | oneBlob.bP.ellipseArea = TMath::Pi()*(rot_wx*rot_wy)/4.0; |
---|
647 | |
---|
648 | // fraction --> (sqrt(2)/2) |
---|
649 | Double_t radiusSquared = 0.; |
---|
650 | if(rot_wx >= rot_wy) |
---|
651 | radiusSquared = (2.*rot_wx*rot_wx) * (0.5/4.); |
---|
652 | else |
---|
653 | radiusSquared = (2.*rot_wy*rot_wy) * (0.5/4.); |
---|
654 | |
---|
655 | oneBlob.bP.circleArea = TMath::Pi()*radiusSquared; |
---|
656 | |
---|
657 | } |
---|
658 | else // no rotation needed |
---|
659 | { |
---|
660 | Double_t radiusSquared = (wx*wx + wy*wy) * 0.5/4.; // fraction --> (sqrt(2)/2) |
---|
661 | oneBlob.bP.circleArea = TMath::Pi()*radiusSquared; |
---|
662 | oneBlob.bP.ellipseArea = TMath::Pi()*(wx * wy)/4.0; |
---|
663 | oneBlob.bP.ellipseA = wx; |
---|
664 | oneBlob.bP.ellipseB = wy; |
---|
665 | } |
---|
666 | |
---|
667 | // finally type initialized to _NOTYPE_ASSIGNED |
---|
668 | // PRBasicSpecies will decide later |
---|
669 | oneBlob.btype = _NOTYPE_ASSIGNED; |
---|
670 | |
---|
671 | //max_x = 0; |
---|
672 | //min_x = GetMatrixXdim()-1; |
---|
673 | //max_y = 0; |
---|
674 | //min_y = GetMatrixYdim()-1; |
---|
675 | return __BLOB_GOOD_TO_STORE; |
---|
676 | } |
---|
677 | |
---|
678 | pair<Int_t, Int_t> BlobsFinder::RotateVector2D(Int_t x, Int_t y, Float_t rotAngle){ |
---|
679 | |
---|
680 | pair<Int_t, Int_t> rotPair; |
---|
681 | |
---|
682 | rotPair.first = |
---|
683 | TMath::FloorNint( (Float_t)x*TMath::Cos(rotAngle) - (Float_t)y*TMath::Sin(rotAngle) ); |
---|
684 | |
---|
685 | rotPair.second = |
---|
686 | TMath::FloorNint( (Float_t)x*TMath::Sin(rotAngle) + (Float_t)y*TMath::Cos(rotAngle) ); |
---|
687 | |
---|
688 | return rotPair; |
---|
689 | } |
---|
690 | |
---|
691 | |
---|
692 | Float_t BlobsFinder::GetBalance(Int_t min_x, Int_t max_x, Int_t min_y, Int_t max_y, |
---|
693 | Float_t * balanceToMin, Float_t * balanceToMax, |
---|
694 | Int_t * minQId, Int_t * maxQId, |
---|
695 | Float_t * slope, Float_t * intersectiony){ |
---|
696 | |
---|
697 | Float_t mean = 0.; |
---|
698 | Int_t divvalmax = 0, divvalmin = 0; |
---|
699 | |
---|
700 | mean = CalcMeanForBlob(min_x, max_x, min_y, max_y, &divvalmax, &divvalmin, minQId, maxQId, slope, intersectiony); |
---|
701 | |
---|
702 | // calculate the width to get the size of one cuadrant, extended to |
---|
703 | // the longest side |
---|
704 | Int_t width_x = (max_x - min_x) + 1; // it is 1 if the width is 1 pixel |
---|
705 | Int_t width_y = (max_y - min_y) + 1; // |
---|
706 | if(width_x < width_y) |
---|
707 | width_x = width_y; |
---|
708 | |
---|
709 | // The balance is the ratio of the minimum ocuppancy to the area |
---|
710 | // of one cuadrant as calculated before. |
---|
711 | *balanceToMin = (Float_t)divvalmin / (Float_t)(width_x*width_x); |
---|
712 | *balanceToMax = (Float_t)divvalmax / mean; |
---|
713 | |
---|
714 | // return mean, probably not used |
---|
715 | return mean; |
---|
716 | } |
---|
717 | |
---|
718 | Float_t BlobsFinder::CalcMeanForBlob(Int_t min_x, Int_t max_x, Int_t min_y, Int_t max_y, |
---|
719 | Int_t * divvalmax, Int_t * divvalmin, |
---|
720 | Int_t * minQId, Int_t * maxQId, |
---|
721 | Float_t * slope, Float_t * intersectiony){ |
---|
722 | |
---|
723 | list< pair < pair<Int_t, Int_t>, Int_t > >::iterator listItr; |
---|
724 | Int_t posx = -1, posy = -1; |
---|
725 | Int_t newposx = -1, newposy = -1; |
---|
726 | Float_t mean = 0.; |
---|
727 | |
---|
728 | ///////////////////////////////////////////////////////////// |
---|
729 | // Working in a doubled grid ! |
---|
730 | min_x *= 2; |
---|
731 | max_x = 2*max_x + 1; |
---|
732 | min_y *= 2; |
---|
733 | max_y = 2*max_y + 1; |
---|
734 | // recalculate centers |
---|
735 | Int_t c_x = TMath::CeilNint((Float_t)(min_x + max_x)/2.0); |
---|
736 | Int_t c_y = TMath::CeilNint((Float_t)(min_y + max_y)/2.0); |
---|
737 | ///////////////////////////////////////////////////////////// |
---|
738 | |
---|
739 | #define NQUADRANTS 4 |
---|
740 | |
---|
741 | Int_t firstQ = 0; |
---|
742 | Int_t secondQ = 0; |
---|
743 | Int_t thirdQ = 0; |
---|
744 | Int_t fourthQ = 0; |
---|
745 | Int_t maxQ = 0; // maximum numbers of counts between the four quadrants |
---|
746 | Int_t minQ = 0; // minimum numbers of counts between the four quadrants |
---|
747 | |
---|
748 | // I need to calculate the mean x and y coordinates in the four cuadrants |
---|
749 | Float_t meansQuadX[] = {0., 0., 0., 0.}; |
---|
750 | Float_t meansQuadY[] = {0., 0., 0., 0.}; |
---|
751 | |
---|
752 | for(listItr = oneBlob.blobContent.begin() ; listItr != oneBlob.blobContent.end() ; listItr++) |
---|
753 | { |
---|
754 | |
---|
755 | posx = (*listItr).first.first; |
---|
756 | posy = (*listItr).first.second; |
---|
757 | |
---|
758 | // Now for each pixel in the list I get 4 |
---|
759 | |
---|
760 | for(Int_t i = 0 ; i < 2 ; i++) |
---|
761 | { |
---|
762 | for(Int_t j = 0 ; j < 2 ; j++) |
---|
763 | { |
---|
764 | // 4 positions in the doubled-grid system |
---|
765 | newposx = 2*posx + i; |
---|
766 | newposy = 2*posy + j; |
---|
767 | |
---|
768 | // now separate quadrants |
---|
769 | if(newposx - c_x < 0) // 3rd or 4th quadrant |
---|
770 | { |
---|
771 | if(newposy - c_y >= 0){ // 4th |
---|
772 | fourthQ++; |
---|
773 | meansQuadX[Q_FOURTH_INDX] += (Float_t)newposx; |
---|
774 | meansQuadY[Q_FOURTH_INDX] += (Float_t)newposy; |
---|
775 | } |
---|
776 | else if(newposy - c_y < 0){ // 3rd |
---|
777 | thirdQ++; |
---|
778 | meansQuadX[Q_THIRD_INDX] += (Float_t)newposx; |
---|
779 | meansQuadY[Q_THIRD_INDX] += (Float_t)newposy; |
---|
780 | } |
---|
781 | } |
---|
782 | else if(newposx - c_x >= 0) // 1st or 2nd quadrant |
---|
783 | { |
---|
784 | if(newposy - c_y >= 0){ // 1st |
---|
785 | firstQ++; |
---|
786 | meansQuadX[Q_FIRST_INDX] += (Float_t)newposx; |
---|
787 | meansQuadY[Q_FIRST_INDX] += (Float_t)newposy; |
---|
788 | } |
---|
789 | else if(newposy - c_y < 0){ // 2nd |
---|
790 | secondQ++; |
---|
791 | meansQuadX[Q_SECOND_INDX] += (Float_t)newposx; |
---|
792 | meansQuadY[Q_SECOND_INDX] += (Float_t)newposy; |
---|
793 | } |
---|
794 | } |
---|
795 | } |
---|
796 | } |
---|
797 | |
---|
798 | } |
---|
799 | |
---|
800 | // means per quadrant |
---|
801 | if(fourthQ) meansQuadX[Q_FOURTH_INDX] /= (Float_t)fourthQ; |
---|
802 | else meansQuadX[Q_FOURTH_INDX] = c_x - 1; // if there is nothing in |
---|
803 | // a cuadrant, take the |
---|
804 | // center +/-1 depending |
---|
805 | // on the cuadrant. This |
---|
806 | // will at least give a |
---|
807 | // correct guess for the |
---|
808 | // slope and ycut |
---|
809 | |
---|
810 | if(fourthQ) meansQuadY[Q_FOURTH_INDX] /= (Float_t)fourthQ; |
---|
811 | else meansQuadY[Q_FOURTH_INDX] = c_y + 1; |
---|
812 | |
---|
813 | if(thirdQ) meansQuadX[Q_THIRD_INDX] /= (Float_t)thirdQ; |
---|
814 | else meansQuadX[Q_THIRD_INDX] = c_x - 1; |
---|
815 | |
---|
816 | if(thirdQ) meansQuadY[Q_THIRD_INDX] /= (Float_t)thirdQ; |
---|
817 | else meansQuadY[Q_THIRD_INDX] = c_y - 1; |
---|
818 | |
---|
819 | if(firstQ) meansQuadX[Q_FIRST_INDX] /= (Float_t)firstQ; |
---|
820 | else meansQuadX[Q_FIRST_INDX] = c_x + 1; |
---|
821 | |
---|
822 | if(firstQ) meansQuadY[Q_FIRST_INDX] /= (Float_t)firstQ; |
---|
823 | else meansQuadY[Q_FIRST_INDX] = c_y + 1; |
---|
824 | |
---|
825 | if(secondQ) meansQuadX[Q_SECOND_INDX] /= (Float_t)secondQ; |
---|
826 | else meansQuadX[Q_SECOND_INDX] = c_x + 1; |
---|
827 | |
---|
828 | if(secondQ) meansQuadY[Q_SECOND_INDX] /= (Float_t)secondQ; |
---|
829 | else meansQuadY[Q_SECOND_INDX] = c_y - 1; |
---|
830 | |
---|
831 | // find max |
---|
832 | if(firstQ > maxQ){ |
---|
833 | maxQ = firstQ; |
---|
834 | *maxQId = Q_FIRST; // Id of the quadrant |
---|
835 | } |
---|
836 | if(secondQ > maxQ){ |
---|
837 | maxQ = secondQ; |
---|
838 | *maxQId = Q_SECOND; |
---|
839 | } |
---|
840 | if(thirdQ > maxQ){ |
---|
841 | maxQ = thirdQ; |
---|
842 | *maxQId = Q_THIRD; |
---|
843 | } |
---|
844 | if(fourthQ > maxQ){ |
---|
845 | maxQ = fourthQ; |
---|
846 | *maxQId = Q_FOURTH; |
---|
847 | } |
---|
848 | // This could be the only option in a perfectly symmetric blob |
---|
849 | minQ = maxQ; |
---|
850 | *minQId = *maxQId; |
---|
851 | |
---|
852 | // find min |
---|
853 | if(firstQ < minQ){ |
---|
854 | minQ = firstQ; |
---|
855 | *minQId = Q_FIRST; |
---|
856 | } |
---|
857 | if(secondQ < minQ){ |
---|
858 | minQ = secondQ; |
---|
859 | *minQId = Q_SECOND; |
---|
860 | } |
---|
861 | if(thirdQ < minQ){ |
---|
862 | minQ = thirdQ; |
---|
863 | *minQId = Q_THIRD; |
---|
864 | } |
---|
865 | if(fourthQ < minQ){ |
---|
866 | minQ = fourthQ; |
---|
867 | *minQId = Q_FOURTH; |
---|
868 | } |
---|
869 | |
---|
870 | // calc mean |
---|
871 | mean = (firstQ + secondQ + thirdQ + fourthQ)/4.0; |
---|
872 | |
---|
873 | // decide the sign of the slope and calculate |
---|
874 | if(*maxQId == Q_FOURTH || *maxQId == Q_SECOND) { // negative slope |
---|
875 | *slope = (meansQuadY[Q_FOURTH_INDX] - meansQuadY[Q_SECOND_INDX])/ |
---|
876 | (meansQuadX[Q_FOURTH_INDX] - meansQuadX[Q_SECOND_INDX]); |
---|
877 | |
---|
878 | *intersectiony = meansQuadY[Q_SECOND_INDX] - ( (*slope) * meansQuadX[Q_SECOND_INDX] ); |
---|
879 | // convert to the simple grid (this working-grid was doubled !) |
---|
880 | *intersectiony /= 2.; |
---|
881 | } |
---|
882 | else if (*maxQId == Q_FIRST || *maxQId == Q_THIRD) { // positive slope |
---|
883 | *slope = (meansQuadY[Q_FIRST_INDX] - meansQuadY[Q_THIRD_INDX])/ |
---|
884 | (meansQuadX[Q_FIRST_INDX] - meansQuadX[Q_THIRD_INDX]); |
---|
885 | |
---|
886 | *intersectiony = meansQuadY[Q_THIRD_INDX] - ( (*slope) * meansQuadX[Q_THIRD_INDX] ); |
---|
887 | // convert to the simple grid (this working-grid was doubled !) |
---|
888 | *intersectiony /= 2.; |
---|
889 | |
---|
890 | } |
---|
891 | |
---|
892 | //Log << MSG::DEBUG << "guess slope = " << *slope << endreq; |
---|
893 | //Log << MSG::DEBUG << "guess ycut = " << *intersectiony << endreq; |
---|
894 | |
---|
895 | // other passed-by-reference values |
---|
896 | *divvalmax = maxQ; |
---|
897 | *divvalmin = minQ; |
---|
898 | |
---|
899 | return mean; |
---|
900 | } |
---|
901 | |
---|
902 | /** |
---|
903 | * Copy properties of all blobs in one frame for ntuple dump |
---|
904 | * |
---|
905 | */ |
---|
906 | void BlobsFinder::CopyPropertiesToNtupleData(){ |
---|
907 | |
---|
908 | if(aB->allBlobs.size() > MAX_BLOBS_PER_FRAME){ |
---|
909 | Log << MSG::ERROR << "Too many blobs found, and can't save to the ntuple. Currently = " << aB->allBlobs.size() << endreq; |
---|
910 | Log << "You need to change MAX_BLOBS_PER_FRAME and recompile this algorithm." << endreq; |
---|
911 | Log << "Or, decide not to save to ntuple at all." << endreq; |
---|
912 | exit(1); |
---|
913 | } |
---|
914 | |
---|
915 | vector<blob>::iterator blobsItr = aB->allBlobs.begin(); |
---|
916 | |
---|
917 | if(m_bP_nBlobs != (Int_t) aB->allBlobs.size()){ // castint uint --> int. Safe. |
---|
918 | Log << MSG::ERROR << "Something went really bad here, check BlobsFinder::CalculateProperties()" << endreq; |
---|
919 | exit(1); |
---|
920 | } |
---|
921 | |
---|
922 | Int_t itr = 0; |
---|
923 | for( ; blobsItr != aB->allBlobs.end() ; blobsItr++) |
---|
924 | { |
---|
925 | m_bP_nPixels[itr] = (*blobsItr).bP.nPixels; |
---|
926 | m_bP_totalCharge[itr] = (*blobsItr).bP.totalCharge; |
---|
927 | |
---|
928 | m_bP_width_x[itr] = (*blobsItr).bP.width_x; |
---|
929 | m_bP_width_y[itr] = (*blobsItr).bP.width_y; |
---|
930 | |
---|
931 | m_bP_geoCenter_x[itr] = (*blobsItr).bP.geoCenter_x; |
---|
932 | m_bP_geoCenter_y[itr] = (*blobsItr).bP.geoCenter_y; |
---|
933 | |
---|
934 | m_bP_weightedCenter_x[itr] = (*blobsItr).bP.weightedCenter_x; |
---|
935 | m_bP_weightedCenter_y[itr] = (*blobsItr).bP.weightedCenter_y; |
---|
936 | |
---|
937 | m_bP_boxArea[itr] = (*blobsItr).bP.boxArea; |
---|
938 | m_bP_circleArea[itr] = (*blobsItr).bP.circleArea; |
---|
939 | m_bP_ellipseArea[itr] = (*blobsItr).bP.ellipseArea; |
---|
940 | |
---|
941 | m_bP_balanceToMin[itr] = (*blobsItr).bP.balanceToMin; |
---|
942 | m_bP_balanceToMax[itr] = (*blobsItr).bP.balanceToMax; |
---|
943 | |
---|
944 | itr++; |
---|
945 | } |
---|
946 | |
---|
947 | } |
---|
948 | |
---|
949 | TString blob::GetTypeAsString(){ |
---|
950 | |
---|
951 | if(btype == _NOTYPE_ASSIGNED) |
---|
952 | return TString("Not_Assigned"); |
---|
953 | else if(btype == _SINGLE_HIT) |
---|
954 | return TString("Single_Hit"); |
---|
955 | else if(btype == _DOUBLE_HIT) |
---|
956 | return TString("Double_Hit"); |
---|
957 | else if(btype == _TRIPLE_HIT) |
---|
958 | return TString("Triple_Hit"); |
---|
959 | else if(btype == _QUAD_HIT) |
---|
960 | return TString("Quad_Hit"); |
---|
961 | else if(btype == _LONG_GAMMA) |
---|
962 | return TString("Long_Gamma"); |
---|
963 | else if(btype == _MIP) |
---|
964 | return TString("Mip"); |
---|
965 | else if(btype == _HEAVY_BLOB) |
---|
966 | return TString("Heavy_Blob"); |
---|
967 | else if(btype == _CURLY) |
---|
968 | return TString("Curly"); |
---|
969 | else if(btype == _UNKNOWN) |
---|
970 | return TString("Unknown"); |
---|
971 | |
---|
972 | return TString("Unknown"); |
---|
973 | } |
---|
974 | |
---|
975 | bool BlobsFinder::SingleHit(){ |
---|
976 | |
---|
977 | if(oneBlob.bP.nPixels == 1) |
---|
978 | return true; |
---|
979 | |
---|
980 | return false; |
---|
981 | } |
---|
982 | |
---|
983 | bool BlobsFinder::DoubleHit(){ |
---|
984 | |
---|
985 | if( |
---|
986 | (oneBlob.bP.nPixels == 2 && oneBlob.bP.width_x == 2 && oneBlob.bP.width_y == 2) |
---|
987 | || |
---|
988 | (oneBlob.bP.nPixels == 2 && oneBlob.bP.width_x == 2 && oneBlob.bP.width_y == 1) |
---|
989 | || |
---|
990 | (oneBlob.bP.nPixels == 2 && oneBlob.bP.width_x == 1 && oneBlob.bP.width_y == 2) |
---|
991 | ) |
---|
992 | return true; |
---|
993 | |
---|
994 | return false; |
---|
995 | } |
---|
996 | |
---|
997 | bool BlobsFinder::TripleHit(){ |
---|
998 | |
---|
999 | if(oneBlob.bP.nPixels == 3 && oneBlob.bP.width_x == 2 && oneBlob.bP.width_y == 2) |
---|
1000 | return true; |
---|
1001 | |
---|
1002 | return false; |
---|
1003 | } |
---|
1004 | |
---|
1005 | bool BlobsFinder::QuadHit(){ |
---|
1006 | |
---|
1007 | if(oneBlob.bP.nPixels == 4 && oneBlob.bP.width_x == 2 && oneBlob.bP.width_y == 2) |
---|
1008 | return true; |
---|
1009 | |
---|
1010 | return false; |
---|
1011 | } |
---|
1012 | |
---|
1013 | void BlobsFinder::CalcMinMaxGeoCenter(Float_t geox, Float_t geoy){ |
---|
1014 | |
---|
1015 | Int_t xcoor = -1; |
---|
1016 | Int_t ycoor = -1; |
---|
1017 | Float_t dist = 0.; |
---|
1018 | // start point |
---|
1019 | xcoor = (*(oneBlob.blobContent.begin())).first.first; |
---|
1020 | ycoor = (*(oneBlob.blobContent.begin())).first.second; |
---|
1021 | dist = TMath::Sqrt((xcoor - geox)*(xcoor - geox) |
---|
1022 | + |
---|
1023 | (ycoor - geoy)*(ycoor - geoy)); |
---|
1024 | oneBlob.bP.minToGeoCenter = dist; |
---|
1025 | oneBlob.bP.maxToGeoCenter = dist; |
---|
1026 | // loop over all elements |
---|
1027 | list< pair < pair<Int_t, Int_t>, Int_t > >::iterator listItr; |
---|
1028 | for(listItr = oneBlob.blobContent.begin() ; listItr != oneBlob.blobContent.end() ; listItr++) |
---|
1029 | { |
---|
1030 | // <x, y> |
---|
1031 | xcoor = (*listItr).first.first; |
---|
1032 | ycoor = (*listItr).first.second; |
---|
1033 | |
---|
1034 | dist = TMath::Sqrt((xcoor - geox)*(xcoor - geox) |
---|
1035 | + |
---|
1036 | (ycoor - geoy)*(ycoor - geoy)); |
---|
1037 | if(dist < oneBlob.bP.minToGeoCenter) |
---|
1038 | oneBlob.bP.minToGeoCenter = dist; |
---|
1039 | else if(dist > oneBlob.bP.maxToGeoCenter) |
---|
1040 | oneBlob.bP.maxToGeoCenter = dist; |
---|
1041 | } |
---|
1042 | |
---|
1043 | } |
---|
1044 | |
---|
1045 | #endif |
---|