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

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