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

Last change on this file since 4555 was 4555, checked in by aliu, 21 years ago
*** empty log message ***
File size: 15.8 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=50;
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 //noise = 0, in the case of MC w/o noise
321 if (noise == 0) noise = 1;
322
323 fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
324
325 }
326
327 //fIsl->SetIslId(fIslId);
328 //fIsl->SetArrivalTime(ftime);
329 fIsl->SetPixNum(fPixNum);
330 fIsl->SetSigToNoise(fSigToNoise);
331 fIsl->SetTimeSpread(timeVariance);
332 fIsl->SetMeanX(meanX);
333 fIsl->SetMeanY(meanY);
334 fIsl->SetDist(dist);
335
336
337 //Length and Width of the larger island according the definition of the hillas parameters
338
339 // calculate 2nd moments
340 // ---------------------
341 Double_t corrxx=0; // [m^2]
342 Double_t corrxy=0; // [m^2]
343 Double_t corryy=0; // [m^2]
344
345 for(Int_t idx=0 ; idx<nPix ; idx++)
346 {
347 MCerPhotPix *pix = fEvt->GetPixById(idx);
348 const MGeomPix &gpix3 = (*fCam)[pix->GetPixId()];
349 const Float_t nphot = pix->GetNumPhotons();
350
351 if (pix == NULL) break;
352
353 if (vect[1][idx]==1){
354
355 const Float_t dx = gpix3.GetX() - X; // [mm]
356 const Float_t dy = gpix3.GetY() - Y; // [mm]
357
358
359 corrxx += nphot * dx*dx; // [mm^2]
360 corrxy += nphot * dx*dy; // [mm^2]
361 corryy += nphot * dy*dy; // [mm^2]
362
363 }
364 }
365
366 // calculate the hillas parameters Width and Length
367
368 const Double_t d0 = corryy - corrxx;
369 const Double_t d1 = corrxy*2;
370 const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1);
371 const Double_t tand = d2 / d1;
372 const Double_t tand2 = tand*tand;
373
374 const Double_t s2 = tand2+1;
375
376 const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/sizeLargeIsl;
377 const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/sizeLargeIsl;
378
379 //
380 // fLength^2 is the second moment along the major axis of the ellipse
381 // fWidth^2 is the second moment along the minor axis of the ellipse
382 //
383 // From the algorithm we get: fWidth <= fLength is always true
384 //
385 // very small numbers can get negative by rounding
386 //
387 length = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm]
388 width = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm]
389
390 fIsl->SetLength(length);
391 fIsl->SetWidth(width);
392
393 // for(Int_t i = 1; i<=numisl ; i++){
394
395 // fIsl->SetDistL(fIsl->GetDist(i-1)/length, i-1);
396 // fIsl->SetDistW(fIsl->GetDist(i-1)/width, i-1);
397 // }
398
399 for(Int_t i = 1; i<=numisl ; i++){
400
401 distL[i-1]=dist[i-1]/length;
402 distW[i-1]=dist[i-1]/width;
403
404 }
405
406 fIsl->SetDistL(distL);
407 fIsl->SetDistW(distW);
408
409 fIsl->SetReadyToSave();
410
411
412 for (Int_t i = 0; i< nVect; i++)
413 delete [] vect[i];
414
415 delete vect;
416
417 return kTRUE;
418}
419
420//------------------------------------------------------------------------------------------
421void MIslandCalc::Calc1(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){
422
423
424 /////////////////////////////
425 //
426 // ALGORITHM # 1
427 // counts the number of islands as you can see in
428 // the event display after doing the std image cleaning
429 //
430 /////////////////////////////
431
432 Int_t sflag;
433 Int_t control;
434
435 Int_t nvect = 0;
436
437 numisl = 0;
438
439 Int_t zeros[nv];
440
441 for(Int_t m = 0; m < nv ; m++)
442 for(Int_t n = 0; n < npix ; n++)
443 vect[m][n] = 0;
444
445 for(Int_t n = 0; n < nv ; n++)
446 zeros[n] = 0;
447
448 //cout << "new event" <<endl;
449 MCerPhotPix *pix;
450
451 //loop over all pixels
452 MCerPhotEvtIter Next(fEvt, kFALSE);
453
454 while ((pix=static_cast<MCerPhotPix*>(Next())))
455 {
456 const Int_t idx = pix->GetPixId();
457
458 const MGeomPix &gpix = (*fCam)[idx];
459 const Int_t nnmax = gpix.GetNumNeighbors();
460
461 if( fEvt->IsPixelUsed(idx))
462 {
463 //cout << idx <<endl;
464 sflag = 0;
465
466 for(Int_t j=0; j < nnmax ; j++)
467 {
468 const Int_t idx2 = gpix.GetNeighbor(j);
469
470 if (idx2 < idx)
471 {
472 for(Int_t k = 1; k <= nvect; k++)
473 {
474 if (vect[k][idx2] == 1)
475 {
476 sflag = 1;
477 vect[k][idx] = 1;
478 }
479 }
480 }
481 }
482
483 if (sflag == 0)
484 {
485 nvect++;
486 // cout << "nvect " << nvect << endl;
487 vect[nvect][idx] = 1;
488 }
489
490 }
491 }
492
493 numisl = nvect;
494
495
496 // Repeated Chain Corrections
497
498
499 for(Int_t i = 1; i <= nvect; i++)
500 {
501 for(Int_t j = i+1; j <= nvect; j++)
502 {
503 control = 0;
504 for(Int_t k = 0; k < npix; k++)
505 {
506 if (vect[i][k] == 1 && vect[j][k] == 1)
507 {
508 control = 1;
509 break;
510 }
511 }
512 if (control == 1)
513 {
514 for(Int_t k = 0; k < npix; k++)
515 {
516 if(vect[j][k] == 1)
517 vect[i][k] = 1;
518 vect[j][k] = 0;
519 zeros[j] = 1;
520 }
521 numisl = numisl-1;
522 }
523
524 }
525 }
526
527
528
529 Int_t l = 1;
530 Int_t numpixels;
531 Int_t pixMAX = 0;
532 Int_t idMAX = 1;
533
534 for(Int_t i = 1; i<= nvect ; i++)
535 {
536 numpixels = 0;
537
538 if (zeros[i] == 0)
539 {
540 for(Int_t k = 0; k<npix; k++)
541 {
542 vect[l][k] = vect[i][k];
543 if (vect[l][k] == 1)
544 numpixels++;
545
546 }
547 if (numpixels>pixMAX)
548 {
549 pixMAX = numpixels;
550 idMAX = l;
551 }
552 l++;
553 }
554 num[i] = numpixels;
555
556 }
557
558 //the larger island will correspond to the 1st component of the vector
559
560 num[nvect +1] = num[1];
561 num[1] = num[idMAX];
562 num[idMAX]=num[1];
563
564 for(Int_t k = 0; k<npix; k++)
565 {
566 vect[nvect+1][k] = vect[1][k];
567 vect[1][k] = vect[idMAX][k];
568 vect[idMAX][k] = vect[nvect+1][k];
569 }
570}
571
572//------------------------------------------------------------------------------------------
573
574void MIslandCalc::Calc2(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){
575
576
577 /////////////////////////////
578 //
579 // ALGORITHM # 2
580 // counts the number of islands considering as the same
581 // islands the ones separated for 2 or less pixels
582 //
583 /////////////////////////////
584
585 Int_t sflag;
586 Int_t control;
587
588 Int_t nvect = 0;
589 numisl = 0;
590
591 Int_t zeros[nv];
592
593 Int_t kk[npix];
594
595 for(Int_t m = 0; m < nv ; m++)
596 for(Int_t n = 0; n < npix ; n++)
597 vect[m][n] = 0;
598
599 for(Int_t n = 0; n < nv ; n++)
600 zeros[n] = 0;
601
602 for(Int_t n = 0; n < npix ; n++)
603 kk[n] = 0;
604
605 MCerPhotPix *pix;
606
607 //1st loop over all pixels
608 MCerPhotEvtIter Next0(fEvt, kFALSE);
609
610 while ((pix=static_cast<MCerPhotPix*>(Next0())))
611 {
612 const Int_t idx = pix->GetPixId();
613
614 const MGeomPix &gpix = (*fCam)[idx];
615 const Int_t nnmax = gpix.GetNumNeighbors();
616
617 if( fEvt->IsPixelUsed(idx))
618 {
619 kk[idx] = 1 ;
620 for(Int_t j=0; j< nnmax; j++)
621 {
622 kk[gpix.GetNeighbor(j)] = 1;
623 }
624 }
625
626 }
627
628 //2nd loop over all pixels
629 MCerPhotEvtIter Next(fEvt, kFALSE);
630
631 while ((pix=static_cast<MCerPhotPix*>(Next())))
632 {
633 const Int_t idx = pix->GetPixId();
634
635 const MGeomPix &gpix = (*fCam)[idx];
636 const Int_t nnmax = gpix.GetNumNeighbors();
637
638 if ( kk[idx] > 0)
639 {
640 sflag = 0;
641
642 for(Int_t j=0; j < nnmax ; j++)
643 {
644 const Int_t idx2 = gpix.GetNeighbor(j);
645
646 if (idx2 < idx)
647 {
648 for(Int_t k = 1; k <= nvect; k++)
649 {
650 if (vect[k][idx2] == 1)
651 {
652 sflag = 1;
653 vect[k][idx] = 1;
654 }
655 }
656 }
657 }
658
659 if (sflag == 0)
660 {
661 nvect++;
662 vect[nvect][idx] = 1;
663 }
664
665 }
666 }
667
668 numisl = nvect;
669
670
671 // Repeated Chain Corrections
672
673 for(Int_t i = 1; i <= nvect; i++)
674 {
675 for(Int_t j = i+1; j <= nvect; j++)
676 {
677 control = 0;
678 for(Int_t k = 0; k < npix; k++)
679 {
680
681 if (vect[i][k] == 1 && vect[j][k] == 1)
682 {
683 control = 1;
684 break;
685 }
686 }
687 if (control == 1)
688 {
689 for(Int_t k = 0; k < npix; k++)
690 {
691 if(vect[j][k] == 1)
692 vect[i][k] = 1;
693 vect[j][k] = 0;
694 zeros[j] = 1;
695 }
696 numisl = numisl-1;
697 }
698
699 }
700 }
701
702
703 Int_t l = 1;
704 Int_t numpixels;
705 Int_t pixMAX = 0;
706 Int_t idMAX = 1;
707
708 for(Int_t i = 1; i<= nvect ; i++)
709 {
710 numpixels = 0;
711
712 if (zeros[i] == 0)
713 {
714 for(Int_t k = 0; k<npix; k++)
715 {
716 vect[l][k] = vect[i][k];
717 if (vect[l][k] == 1)
718 numpixels++;
719 }
720 if (numpixels>pixMAX)
721 {
722 pixMAX = numpixels;
723 idMAX = l;
724 }
725 l++;
726 }
727 num[i] = numpixels;
728 }
729
730
731 //the larger island will correspond to the 1st component of the vector
732
733 num[nvect +1] = num[1];
734 num[1] = num[idMAX];
735 num[idMAX]=num[1];
736
737 for(Int_t k = 0; k<npix; k++)
738 {
739 vect[nvect+1][k] = vect[1][k];
740 vect[1][k] = vect[idMAX][k];
741 vect[idMAX][k] = vect[nvect+1][k];
742 }
743
744}
745
Note: See TracBrowser for help on using the repository browser.