source: trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc@ 5258

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