source: trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandCalc.cc@ 4450

Last change on this file since 4450 was 4442, checked in by aliu, 20 years ago
*** empty log message ***
File size: 15.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Ester Aliu, 2/2004 <aliu@ifae.es>
19|
20! Last Update: 7/2004
21!
22!
23! Copyright: MAGIC Software Development, 2000-2004
24!
25!
26\* ======================================================================== */
27
28/////////////////////////////////////////////////////////////////////////////
29//
30// MIslandCalc
31//
32// The Island Calc task calculates the some islands parameters for each
33// of the events such as:
34//
35// - fIslNum // number of islands
36// - fIslId[577] // island Id
37// - fPixNum[numisl] // number of pixels in the island
38// - fSigToNoise[numisl] // signal to noise of the island
39// - ftime[numisl][577] // mean of the arrival time
40// - fTimeSpread[numisl] // mean arrival time spread of the island
41// - fMeanX[numisl] // mean X position of the island
42// - fMeanY[numisl] // mean Y position of the island
43// - fDist[numisl] // dist between an island and the continent
44// - fLength // major axis of the larger island ellipse
45// - fWidth // minor axis of the larger island ellipse
46// - fDistL[numisl] // dist divided by lenght of the larger island
47// - fDistW[numisl] // dist divided by width of the larger island
48//
49// Input Containers:
50// MGeomCam
51// MCerPhotEvt
52// MPedestalCam
53// MArrivalTimeCam
54//
55// Output Containers:
56// MIslands
57//
58/////////////////////////////////////////////////////////////////////////////
59#include "MIslandCalc.h"
60
61#include <stdlib.h> // atof
62#include <fstream> // ofstream, SavePrimitive
63
64#include "MLog.h"
65#include "MLogManip.h"
66
67#include "MIslands.h"
68
69#include "MParList.h"
70
71#include "MGeomPix.h"
72#include "MGeomCam.h"
73
74#include "MCerPhotPix.h"
75#include "MCerPhotEvt.h"
76
77#include "MPedestalCam.h"
78#include "MPedestalPix.h"
79
80#include "MArrivalTimeCam.h"
81#include "MArrivalTimePix.h"
82
83ClassImp(MIslandCalc);
84
85
86using namespace std;
87
88// --------------------------------------------------------------------------
89//
90// Default constructor.
91//
92MIslandCalc::MIslandCalc(const char* name, const char* title)
93 : fIsl(NULL)
94{
95 fName = name ? name : "MIslandCalc";
96 fTitle = title ? title : "Calculate island parameters";
97}
98
99
100// --------------------------------------------------------------------------
101Int_t MIslandCalc::PreProcess (MParList *pList)
102{
103 fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
104 if (!fCam)
105 {
106 *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
107 return kFALSE;
108 }
109
110 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
111 if (!fEvt)
112 {
113 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
114 return kFALSE;
115 }
116
117 fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
118 if (!fPed)
119 {
120 *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
121 return kFALSE;
122 }
123
124 fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam"));
125 if (!fTime)
126 {
127 *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl;
128 return kFALSE;
129 }
130
131 if (strlen(fIslName) > 0)
132 {
133 fIsl = (MIslands*)pList->FindCreateObj("MIslands", AddSerialNumber(fIslName));
134 //cout << "kk1" << endl;
135 }
136 else
137 {
138 fIsl = (MIslands*)pList->FindCreateObj(AddSerialNumber("MIslands"));
139 //cout << "kk2" << endl;
140 }
141 if (!fIsl)
142 return kFALSE;
143
144 return kTRUE;
145}
146
147
148Int_t MIslandCalc::Process(){
149
150 IslandPar();
151
152 return kTRUE;
153
154}
155
156
157Int_t MIslandCalc::IslandPar(){
158
159 //calculates all the island parameters
160
161 const Int_t nPix=577;
162 const Int_t nVect=20;
163 Int_t numisl;
164
165 Int_t** vect;
166 vect = new Int_t*[nVect];
167 for(Int_t i=0;i<nVect;i++)
168 vect[i]= new Int_t[nPix];
169
170 Int_t num[nVect];
171 // num = new Int_t*[nVect];
172
173 if (fIslandAlgorithm == 1)
174 Calc1(numisl,nVect,nPix,vect,num);
175 if (fIslandAlgorithm == 2)
176 Calc2(numisl,nVect,nPix,vect,num);
177
178 //set the number of islands in one event
179
180 fIsl->SetIslNum(numisl);
181
182
183 //examine each island...
184
185 Float_t noise;
186 Float_t signal;
187
188 /*
189 Float_t** ftime;
190 ftime = new Float_t*[numisl];
191
192 Int_t pixNumIsl = 0;
193
194 for(Int_t i=0;i<numisl;i++)
195 {
196 pixNumIsl = num[i+1];
197 ftime[i]= new Float_t[pixNumIsl];
198 }
199
200 Int_t** fIslId;
201 fIslId = new Int_t*[numisl];
202 for(Int_t i=0;i<numisl;i++)
203 {
204 pixNumIsl = num[i+1];
205 fIslId[i]= new Int_t[pixNumIsl];
206 }
207 */
208
209 Int_t fPixNum[numisl];
210 Float_t fSigToNoise[numisl];
211 Float_t time[nPix];
212 Float_t timeVariance[numisl];
213 Float_t meanX[numisl];
214 Float_t meanY[numisl];
215 Float_t dist[numisl];
216 Float_t distL[numisl];
217 Float_t distW[numisl];
218 Float_t size, sizeLargeIsl, length, width;
219
220
221 //reset the "sets" functions
222 if (numisl <1)
223 fIsl->SetIslNum(0);
224
225 for(Int_t i = 0; i<10 ;i++){
226 for(Int_t idx = 0; idx<nPix; idx++)
227 {
228 if (i == 0)
229 fIsl->SetIslId(idx, -1);
230 //fIsl->SetArrivalTime(i, idx, -1 );
231 }
232 }
233
234 Float_t X = 0;
235 Float_t Y = 0;
236 sizeLargeIsl = 0;
237
238
239 for(Int_t i = 1; i<=numisl ; i++)
240 {
241 Int_t n = 0;
242 //Int_t ncore = 0;
243
244 Float_t MIN = 10000.;
245 Float_t MAX = 0.;
246
247 signal = 0;
248 noise = 0;
249
250 size = 0;
251 meanX[i-1] = 0;
252 meanY[i-1] = 0;
253 dist[i-1] = 0;
254
255 fPixNum[i-1] = 0;
256 timeVariance[i-1] = 0;
257
258
259 for(Int_t idx=0 ; idx<nPix ; idx++)
260 {
261
262 MCerPhotPix *pix = fEvt->GetPixById(idx);
263 const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()];
264 const MPedestalPix &ped = (*fPed)[idx];
265 const MArrivalTimePix &timepix = (*fTime)[idx];
266 const Float_t nphot = pix->GetNumPhotons();
267
268 if (pix == NULL) break;
269
270 if (vect[i][idx]==1){
271
272 fPixNum[i-1]++;
273 signal += nphot * (fCam->GetPixRatio(idx));
274 noise += pow(ped.GetPedestalRms(),2);
275
276 size += nphot;
277 if (i == 1)
278 sizeLargeIsl += nphot;
279
280 meanX[i-1] += nphot * gpix2.GetX();
281 meanY[i-1] += nphot * gpix2.GetY();
282
283 time[i-1] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
284
285 // ftime[i-1][n] = time[i-1];
286 // fIslId[i-1][n] = idx;
287
288 //calculates the time spread only for core pixels
289 if (fEvt->IsPixelCore(idx)){
290
291 if (time[i-1] > MAX)
292 MAX = time[i-1];
293 if (time[i-1] < MIN)
294 MIN = time[i-1];
295 // ncore++;
296 }
297
298 fIsl->SetIslId(idx, i-1);
299 // fIsl->SetArrivalTime(i-1, idx, time[n]);
300
301 n++;
302 }
303
304 }
305
306 meanX[i-1] /= size;
307 meanY[i-1] /= size;
308
309
310 if (i == 1){
311 X = meanX[i-1];
312 Y = meanY[i-1];
313 }
314
315 dist[i-1] = TMath::Power(meanX[i-1]-X,2) + TMath::Power(meanY[i-1]-Y,2);
316 dist[i-1] = TMath::Sqrt(dist[i-1]);
317
318 timeVariance[i-1] = MAX-MIN;
319
320 fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
321
322 }
323
324 //fIsl->SetIslId(fIslId);
325 //fIsl->SetArrivalTime(ftime);
326 fIsl->SetPixNum(fPixNum);
327 fIsl->SetSigToNoise(fSigToNoise);
328 fIsl->SetTimeSpread(timeVariance);
329 fIsl->SetMeanX(meanX);
330 fIsl->SetMeanY(meanY);
331 fIsl->SetDist(dist);
332
333
334 //Length and Width of the larger island according the definition of the hillas parameters
335
336 // calculate 2nd moments
337 // ---------------------
338 Double_t corrxx=0; // [m^2]
339 Double_t corrxy=0; // [m^2]
340 Double_t corryy=0; // [m^2]
341
342 for(Int_t idx=0 ; idx<nPix ; idx++)
343 {
344 MCerPhotPix *pix = fEvt->GetPixById(idx);
345 const MGeomPix &gpix3 = (*fCam)[pix->GetPixId()];
346 const Float_t nphot = pix->GetNumPhotons();
347
348 if (pix == NULL) break;
349
350 if (vect[1][idx]==1){
351
352 const Float_t dx = gpix3.GetX() - X; // [mm]
353 const Float_t dy = gpix3.GetY() - Y; // [mm]
354
355
356 corrxx += nphot * dx*dx; // [mm^2]
357 corrxy += nphot * dx*dy; // [mm^2]
358 corryy += nphot * dy*dy; // [mm^2]
359
360 }
361 }
362
363 // calculate the hillas parameters Width and Length
364
365 const Double_t d0 = corryy - corrxx;
366 const Double_t d1 = corrxy*2;
367 const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1);
368 const Double_t tand = d2 / d1;
369 const Double_t tand2 = tand*tand;
370
371 const Double_t s2 = tand2+1;
372
373 const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/sizeLargeIsl;
374 const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/sizeLargeIsl;
375
376 //
377 // fLength^2 is the second moment along the major axis of the ellipse
378 // fWidth^2 is the second moment along the minor axis of the ellipse
379 //
380 // From the algorithm we get: fWidth <= fLength is always true
381 //
382 // very small numbers can get negative by rounding
383 //
384 length = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm]
385 width = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm]
386
387 fIsl->SetLength(length);
388 fIsl->SetWidth(width);
389
390 // for(Int_t i = 1; i<=numisl ; i++){
391
392 // fIsl->SetDistL(fIsl->GetDist(i-1)/length, i-1);
393 // fIsl->SetDistW(fIsl->GetDist(i-1)/width, i-1);
394 // }
395
396 for(Int_t i = 1; i<=numisl ; i++){
397
398 distL[i-1]=dist[i-1]/length;
399 distW[i-1]=dist[i-1]/width;
400
401 }
402
403 fIsl->SetDistL(distL);
404 fIsl->SetDistW(distW);
405
406 fIsl->SetReadyToSave();
407
408
409 for(Int_t i=0;i<nVect;i++)
410 delete [] vect[i];
411
412 delete vect;
413
414 return kTRUE;
415}
416
417//------------------------------------------------------------------------------------------
418void MIslandCalc::Calc1(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){
419
420
421 /////////////////////////////
422 //
423 // ALGORITHM # 1
424 // counts the number of islands as you can see in
425 // the event display after doing the std image cleaning
426 //
427 /////////////////////////////
428
429 Int_t sflag;
430 Int_t control;
431
432 Int_t nvect = 0;
433
434 numisl = 0;
435
436 Int_t zeros[nv];
437
438 for(Int_t m = 0; m < nv ; m++)
439 for(Int_t n = 0; n < npix ; n++)
440 vect[m][n] = 0;
441
442 for(Int_t n = 0; n < nv ; n++)
443 zeros[n] = 0;
444
445
446 MCerPhotPix *pix;
447
448 //loop over all pixels
449 MCerPhotEvtIter Next(fEvt, kFALSE);
450
451 while ((pix=static_cast<MCerPhotPix*>(Next())))
452 {
453 const Int_t idx = pix->GetPixId();
454
455 const MGeomPix &gpix = (*fCam)[idx];
456 const Int_t nnmax = gpix.GetNumNeighbors();
457
458 if( fEvt->IsPixelUsed(idx))
459 {
460 sflag = 0;
461
462 for(Int_t j=0; j < nnmax ; j++)
463 {
464 const Int_t idx2 = gpix.GetNeighbor(j);
465
466 if (idx2 < idx)
467 {
468 for(Int_t k = 1; k <= nvect; k++)
469 {
470 if (vect[k][idx2] == 1)
471 {
472 sflag = 1;
473 vect[k][idx] = 1;
474 }
475 }
476 }
477 }
478
479 if (sflag == 0)
480 {
481 nvect++;
482 vect[nvect][idx] = 1;
483 }
484
485 }
486 }
487
488 numisl = nvect;
489
490
491 // Repeated Chain Corrections
492
493
494 for(Int_t i = 1; i <= nvect; i++)
495 {
496 for(Int_t j = i+1; j <= nvect; j++)
497 {
498 control = 0;
499 for(Int_t k = 0; k < npix; k++)
500 {
501 if (vect[i][k] == 1 && vect[j][k] == 1)
502 {
503 control = 1;
504 break;
505 }
506 }
507 if (control == 1)
508 {
509 for(Int_t k = 0; k < npix; k++)
510 {
511 if(vect[j][k] == 1)
512 vect[i][k] = 1;
513 vect[j][k] = 0;
514 zeros[j] = 1;
515 }
516 numisl = numisl-1;
517 }
518
519 }
520 }
521
522
523
524 Int_t l = 1;
525 Int_t numpixels;
526 Int_t pixMAX = 0;
527 Int_t idMAX = 1;
528
529 for(Int_t i = 1; i<= nvect ; i++)
530 {
531 numpixels = 0;
532
533 if (zeros[i] == 0)
534 {
535 for(Int_t k = 0; k<npix; k++)
536 {
537 vect[l][k] = vect[i][k];
538 if (vect[l][k] == 1)
539 numpixels++;
540
541 }
542 if (numpixels>pixMAX)
543 {
544 pixMAX = numpixels;
545 idMAX = l;
546 }
547 l++;
548 }
549 num[i] = numpixels;
550
551 }
552
553 //the larger island will correspond to the 1st component of the vector
554
555 num[nvect +1] = num[1];
556 num[1] = num[idMAX];
557 num[idMAX]=num[1];
558
559 for(Int_t k = 0; k<npix; k++)
560 {
561 vect[nvect+1][k] = vect[1][k];
562 vect[1][k] = vect[idMAX][k];
563 vect[idMAX][k] = vect[nvect+1][k];
564 }
565}
566
567//------------------------------------------------------------------------------------------
568
569void MIslandCalc::Calc2(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){
570
571
572 /////////////////////////////
573 //
574 // ALGORITHM # 2
575 // counts the number of islands considering as the same
576 // islands the ones separated for 2 or less pixels
577 //
578 /////////////////////////////
579
580 Int_t sflag;
581 Int_t control;
582
583 Int_t nvect = 0;
584 numisl = 0;
585
586 Int_t zeros[nv];
587
588 Int_t kk[npix];
589
590 for(Int_t m = 0; m < nv ; m++)
591 for(Int_t n = 0; n < npix ; n++)
592 vect[m][n] = 0;
593
594 for(Int_t n = 0; n < nv ; n++)
595 zeros[n] = 0;
596
597 for(Int_t n = 0; n < npix ; n++)
598 kk[n] = 0;
599
600 MCerPhotPix *pix;
601
602 //1st loop over all pixels
603 MCerPhotEvtIter Next0(fEvt, kFALSE);
604
605 while ((pix=static_cast<MCerPhotPix*>(Next0())))
606 {
607 const Int_t idx = pix->GetPixId();
608
609 const MGeomPix &gpix = (*fCam)[idx];
610 const Int_t nnmax = gpix.GetNumNeighbors();
611
612 if( fEvt->IsPixelUsed(idx))
613 {
614 kk[idx] = 1 ;
615 for(Int_t j=0; j< nnmax; j++)
616 {
617 kk[gpix.GetNeighbor(j)] = 1;
618 }
619 }
620
621 }
622
623 //2nd loop over all pixels
624 MCerPhotEvtIter Next(fEvt, kFALSE);
625
626 while ((pix=static_cast<MCerPhotPix*>(Next())))
627 {
628 const Int_t idx = pix->GetPixId();
629
630 const MGeomPix &gpix = (*fCam)[idx];
631 const Int_t nnmax = gpix.GetNumNeighbors();
632
633 if ( kk[idx] > 0)
634 {
635 sflag = 0;
636
637 for(Int_t j=0; j < nnmax ; j++)
638 {
639 const Int_t idx2 = gpix.GetNeighbor(j);
640
641 if (idx2 < idx)
642 {
643 for(Int_t k = 1; k <= nvect; k++)
644 {
645 if (vect[k][idx2] == 1)
646 {
647 sflag = 1;
648 vect[k][idx] = 1;
649 }
650 }
651 }
652 }
653
654 if (sflag == 0)
655 {
656 nvect++;
657 vect[nvect][idx] = 1;
658 }
659
660 }
661 }
662
663 numisl = nvect;
664
665
666 // Repeated Chain Corrections
667
668 for(Int_t i = 1; i <= nvect; i++)
669 {
670 for(Int_t j = i+1; j <= nvect; j++)
671 {
672 control = 0;
673 for(Int_t k = 0; k < npix; k++)
674 {
675
676 if (vect[i][k] == 1 && vect[j][k] == 1)
677 {
678 control = 1;
679 break;
680 }
681 }
682 if (control == 1)
683 {
684 for(Int_t k = 0; k < npix; k++)
685 {
686 if(vect[j][k] == 1)
687 vect[i][k] = 1;
688 vect[j][k] = 0;
689 zeros[j] = 1;
690 }
691 numisl = numisl-1;
692 }
693
694 }
695 }
696
697
698 Int_t l = 1;
699 Int_t numpixels;
700 Int_t pixMAX = 0;
701 Int_t idMAX = 1;
702
703 for(Int_t i = 1; i<= nvect ; i++)
704 {
705 numpixels = 0;
706
707 if (zeros[i] == 0)
708 {
709 for(Int_t k = 0; k<npix; k++)
710 {
711 vect[l][k] = vect[i][k];
712 if (vect[l][k] == 1)
713 numpixels++;
714 }
715 if (numpixels>pixMAX)
716 {
717 pixMAX = numpixels;
718 idMAX = l;
719 }
720 l++;
721 }
722 num[i] = numpixels;
723 }
724
725
726 //the larger island will correspond to the 1st component of the vector
727
728 num[nvect +1] = num[1];
729 num[1] = num[idMAX];
730 num[idMAX]=num[1];
731
732 for(Int_t k = 0; k<npix; k++)
733 {
734 vect[nvect+1][k] = vect[1][k];
735 vect[1][k] = vect[idMAX][k];
736 vect[idMAX][k] = vect[nvect+1][k];
737 }
738
739}
740
Note: See TracBrowser for help on using the repository browser.