source: Idarraga/mafalda/BlobsFinder/BlobsFinder.cpp @ 163

Last change on this file since 163 was 163, checked in by idarraga, 13 years ago
File size: 28.4 KB
Line 
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
14using namespace MSG;
15
16ClassImp(BlobsFinder)
17ClassImp(AllBlobsContainer)
18
19/*
20 * AllBlobsContainer constructor
21 *
22 */
23AllBlobsContainer::AllBlobsContainer(MediPixAlgo * algo) : 
24CandidateContainer(algo) {
25
26}
27
28void AllBlobsContainer::SetBlobType(Int_t id, blobtype type){
29        allBlobs[id].SetBlobType(type);
30}
31
32blobtype AllBlobsContainer::GetBlobType(Int_t id) { 
33        return allBlobs[id].GetBlobType();
34}
35
36TString AllBlobsContainer::GetTypeAsString(Int_t id) { 
37        return allBlobs[id].GetTypeAsString();
38}
39
40void blob::SetBlobType(blobtype type)
41{
42        btype = type;
43}
44
45/**
46 * BlobsFinder constructor
47 *
48 */
49BlobsFinder::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 */
79Int_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 */
93void 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 */
112void 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 */
146void 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 */
223void 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 */
232void 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 */
259bool 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 */
292bool 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 */
352bool 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 */
366void 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
381void 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 */
398void 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 */
412pair<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 */
449void 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 */
465int 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
678pair<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
692Float_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
718Float_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 */
906void 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
949TString 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
975bool BlobsFinder::SingleHit(){
976
977        if(oneBlob.bP.nPixels == 1)
978                return true;
979
980        return false;
981}
982
983bool 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
997bool 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
1005bool 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
1013void 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
Note: See TracBrowser for help on using the repository browser.