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

Last change on this file since 4406 was 4286, checked in by aliu, 21 years ago
*** empty log message ***
File size: 11.6 KB
Line 
1/* ======================================================================== *\
2!
3!
4! Author(s): Ester Aliu, 2/2004 <aliu@ifae.es>
5!
6! Last Update: 6/2004
7!
8\* ======================================================================== */
9
10/////////////////////////////////////////////////////////////////////////////
11//
12// MIslandCalc
13//
14// The Island Calc task calculates the some islands parameters for each
15// of the events such as:
16//
17// - fIslNum // number of islands
18// - fIslId[577] // island Id
19// - fPixNum[20]; // number of pixels in the island
20// - fSigToNoise[20]; // signal to noise of the island
21// - fTime[20][577] // mean of the arrival time
22// - fTimeSpread[20]; // mean arrival time spread of the island
23//
24//
25// Input Containers:
26// MGeomCam
27// MCerPhotEvt
28// MPedestalCam
29// MArrivalTimeCam
30//
31// Output Containers:
32// MIslands
33//
34/////////////////////////////////////////////////////////////////////////////
35#include "MIslandCalc.h"
36
37#include <stdlib.h> // atof
38#include <fstream> // ofstream, SavePrimitive
39
40#include "MLog.h"
41#include "MLogManip.h"
42
43#include "MIslands.h"
44
45#include "MParList.h"
46
47#include "MGeomPix.h"
48#include "MGeomCam.h"
49
50#include "MCerPhotPix.h"
51#include "MCerPhotEvt.h"
52
53#include "MPedestalCam.h"
54#include "MPedestalPix.h"
55
56#include "MArrivalTimeCam.h"
57#include "MArrivalTimePix.h"
58
59ClassImp(MIslandCalc);
60
61
62using namespace std;
63
64// --------------------------------------------------------------------------
65//
66// Default constructor.
67//
68MIslandCalc::MIslandCalc(const char* name, const char* title)
69 : fIsl(NULL)
70{
71 fName = name ? name : "MIslandCalc";
72 fTitle = title ? title : "Calculate island parameters";
73}
74
75
76// --------------------------------------------------------------------------
77Int_t MIslandCalc::PreProcess (MParList *pList)
78{
79 fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
80 if (!fCam)
81 {
82 *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
83 return kFALSE;
84 }
85
86 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
87 if (!fEvt)
88 {
89 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
90 return kFALSE;
91 }
92
93 fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
94 if (!fPed)
95 {
96 *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
97 return kFALSE;
98 }
99
100 fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam"));
101 if (!fTime)
102 {
103 *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl;
104 return kFALSE;
105 }
106
107 if (strlen(fIslName) > 0)
108 {
109 fIsl = (MIslands*)pList->FindCreateObj("MIslands", AddSerialNumber(fIslName));
110 //cout << "kk1" << endl;
111 }
112 else
113 {
114 fIsl = (MIslands*)pList->FindCreateObj(AddSerialNumber("MIslands"));
115 //cout << "kk2" << endl;
116 }
117 if (!fIsl)
118 return kFALSE;
119
120 return kTRUE;
121}
122
123
124Int_t MIslandCalc::Process()
125{
126
127 if (fIslandAlgorithm == 1)
128 Calc1();
129
130 if (fIslandAlgorithm == 2)
131 Calc2();
132 return kTRUE;
133
134}
135
136
137Int_t MIslandCalc::Calc1(){
138
139
140 /////////////////////////////
141 //
142 // ALGORITHM # 1
143 // counts the number of algorithms as you can see in
144 // the event display after doing the std image cleaning
145 //
146 /////////////////////////////
147
148 Float_t noise;
149 Float_t signal;
150
151 Int_t npix = 577;
152
153 Int_t sflag;
154 Int_t control;
155
156 Int_t nvect = 0;
157 Int_t fIslNum = 0;
158
159 Int_t vect[50][577];
160
161 Int_t zeros[50];
162
163 for(Int_t m = 0; m < 50 ; m++)
164 for(Int_t n = 0; n < npix ; n++)
165 vect[m][n] = 0;
166
167 for(Int_t n = 0; n < 50 ; n++)
168 zeros[n] = 0;
169
170
171 MCerPhotPix *pix;
172
173 //loop over all pixels
174 MCerPhotEvtIter Next(fEvt, kFALSE);
175
176 while ((pix=static_cast<MCerPhotPix*>(Next())))
177 {
178 const Int_t idx = pix->GetPixId();
179
180 const MGeomPix &gpix = (*fCam)[idx];
181 const Int_t nnmax = gpix.GetNumNeighbors();
182
183 if( fEvt->IsPixelUsed(idx))
184 {
185 sflag = 0;
186
187 for(Int_t j=0; j < nnmax ; j++)
188 {
189 const Int_t idx2 = gpix.GetNeighbor(j);
190
191 if (idx2 < idx)
192 {
193 for(Int_t k = 1; k <= nvect; k++)
194 {
195 if (vect[k][idx2] == 1)
196 {
197 sflag = 1;
198 vect[k][idx] = 1;
199 }
200 }
201 }
202 }
203
204 if (sflag == 0)
205 {
206 nvect++;
207 vect[nvect][idx] = 1;
208 }
209
210 }
211 }
212
213 fIslNum = nvect;
214
215
216 // Repeated Chain Corrections
217
218 for(Int_t i = 1; i <= nvect; i++)
219 {
220 for(Int_t j = i+1; j <= nvect; j++)
221 {
222 control = 0;
223 for(Int_t k = 0; k < npix; k++)
224 {
225
226 if (vect[i][k] == 1 && vect[j][k] == 1)
227 {
228 control = 1;
229 break;
230 }
231 }
232 if (control == 1)
233 {
234 for(Int_t k = 0; k < npix; k++)
235 {
236 if(vect[j][k] == 1)
237 vect[i][k] = 1;
238 vect[j][k] = 0;
239 zeros[j] = 1;
240 }
241 fIslNum = fIslNum-1;
242 }
243
244 }
245 }
246
247
248 Int_t l = 1;
249
250 for(Int_t i = 1; i<= nvect ; i++)
251 {
252 if (zeros[i] == 0)
253 {
254 for(Int_t k = 0; k<npix; k++)
255 {
256 vect[l][k] = vect[i][k];
257 }
258 l++;
259 }
260 }
261
262
263 //set the number of islands in one event
264 fIsl->SetIslNum(fIslNum);
265
266 //examine each island...
267 Int_t fPixNum[fIslNum];
268 Float_t fSigToNoise[fIslNum];
269 Float_t time[577];
270 Float_t timeVariance[fIslNum];
271
272 //reset the "sets" functions
273 if (fIslNum <1)
274 fIsl->SetIslNum(0);
275
276 for(Int_t i = 0; i<20 ;i++)
277 {
278 fIsl->SetPixNum(i,-1);
279 fIsl->SetSigToNoise(i,-1);
280 fIsl->SetTimeSpread(i,-1);
281
282 for(Int_t idx = 0; idx<npix; idx++)
283 {
284 fIsl->SetIslId(idx, -1);
285 fIsl->SetArrivalTime(i, idx, -1 );
286 }
287 }
288
289
290 for(Int_t i = 1; i<=fIslNum ; i++)
291 {
292 Int_t n = 0;
293 Int_t ncore = 0;
294
295 Float_t MIN = 10000;
296 Float_t MAX = 0;
297
298 signal = 0;
299 noise = 0;
300 fPixNum[i-1] = 0;
301 timeVariance[i-1] = 0;
302
303 for(Int_t idx=0 ; idx<npix ; idx++)
304 {
305
306 MCerPhotPix *pix = fEvt->GetPixById(idx);
307 const MPedestalPix &ped = (*fPed)[idx];
308 const MArrivalTimePix &timepix = (*fTime)[idx];
309
310 if (pix == NULL) break;
311
312 if (vect[i][idx]==1){
313
314 fPixNum[i-1]++;
315 signal += pix->GetNumPhotons() * (fCam->GetPixRatio(idx));
316 noise += pow(ped.GetPedestalRms(),2);
317
318 time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
319
320
321 if (fEvt->IsPixelCore(idx)){
322
323 if (time[n] > MAX)
324 MAX = time[n];
325 if (time[n] < MIN)
326 MIN = time[n];
327
328 ncore++;
329 }
330
331 fIsl->SetIslId(idx, i-1);
332 fIsl->SetArrivalTime(i-1, idx, time[n]);
333
334 n++;
335 }
336
337 }
338
339 // Float_t mean = timeVariance[i-1]/ncore;
340
341 timeVariance[i-1] = 0;
342
343 timeVariance[i-1] = (MAX - MIN)/ncore;
344 timeVariance[i-1] = (MAX - MIN)/ncore;
345
346 fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
347
348 fIsl->SetPixNum(i-1,fPixNum[i-1]);
349 fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
350 fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
351
352 }
353
354 fIsl->SetReadyToSave();
355
356 return kTRUE;
357}
358
359
360Int_t MIslandCalc::Calc2(){
361
362
363 /////////////////////////////
364 //
365 // ALGORITHM # 2
366 // counts the number of islands considering as the same
367 // islands the ones separated for 2 or less pixels
368 //
369 /////////////////////////////
370
371 Float_t noise;
372 Float_t signal;
373
374 Int_t npix = 577;
375
376 Int_t sflag;
377 Int_t control;
378
379 Int_t nvect = 0;
380 Int_t fIslNum = 0;
381
382 Int_t vect[50][577];
383
384 Int_t zeros[50];
385
386 Int_t kk[577];
387
388 for(Int_t m = 0; m < 50 ; m++)
389 for(Int_t n = 0; n < npix ; n++)
390 vect[m][n] = 0;
391
392 for(Int_t n = 0; n < 50 ; n++)
393 zeros[n] = 0;
394
395 for(Int_t n = 0; n < 577 ; n++)
396 kk[n] = 0;
397
398 MCerPhotPix *pix;
399
400 //1st loop over all pixels
401 MCerPhotEvtIter Next0(fEvt, kFALSE);
402
403 while ((pix=static_cast<MCerPhotPix*>(Next0())))
404 {
405 const Int_t idx = pix->GetPixId();
406
407 const MGeomPix &gpix = (*fCam)[idx];
408 const Int_t nnmax = gpix.GetNumNeighbors();
409
410 if( fEvt->IsPixelUsed(idx))
411 {
412 kk[idx] = 1 ;
413 for(Int_t j=0; j< nnmax; j++)
414 {
415 kk[gpix.GetNeighbor(j)] = 1;
416 }
417 }
418
419 }
420
421 //2nd loop over all pixels
422 MCerPhotEvtIter Next(fEvt, kFALSE);
423
424 while ((pix=static_cast<MCerPhotPix*>(Next())))
425 {
426 const Int_t idx = pix->GetPixId();
427
428 const MGeomPix &gpix = (*fCam)[idx];
429 const Int_t nnmax = gpix.GetNumNeighbors();
430
431 if ( kk[idx] > 0)
432 {
433 sflag = 0;
434
435 for(Int_t j=0; j < nnmax ; j++)
436 {
437 const Int_t idx2 = gpix.GetNeighbor(j);
438
439 if (idx2 < idx)
440 {
441 for(Int_t k = 1; k <= nvect; k++)
442 {
443 if (vect[k][idx2] == 1)
444 {
445 sflag = 1;
446 vect[k][idx] = 1;
447 }
448 }
449 }
450 }
451
452 if (sflag == 0)
453 {
454 nvect++;
455 vect[nvect][idx] = 1;
456 }
457
458 }
459 }
460
461 fIslNum = nvect;
462
463
464 // Repeated Chain Corrections
465
466 for(Int_t i = 1; i <= nvect; i++)
467 {
468 for(Int_t j = i+1; j <= nvect; j++)
469 {
470 control = 0;
471 for(Int_t k = 0; k < npix; k++)
472 {
473
474 if (vect[i][k] == 1 && vect[j][k] == 1)
475 {
476 control = 1;
477 break;
478 }
479 }
480 if (control == 1)
481 {
482 for(Int_t k = 0; k < npix; k++)
483 {
484 if(vect[j][k] == 1)
485 vect[i][k] = 1;
486 vect[j][k] = 0;
487 zeros[j] = 1;
488 }
489 fIslNum = fIslNum-1;
490 }
491
492 }
493 }
494
495
496 Int_t l = 1;
497
498 for(Int_t i = 1; i<= nvect ; i++)
499 {
500 if (zeros[i] == 0)
501 {
502 for(Int_t k = 0; k<npix; k++)
503 {
504 vect[l][k] = vect[i][k];
505 }
506 l++;
507 }
508 }
509
510
511 //set the number of islands in one event
512 fIsl->SetIslNum(fIslNum);
513
514 //examine each island...
515 Int_t fPixNum[fIslNum];
516 Float_t fSigToNoise[fIslNum];
517 Float_t time[577];
518 Float_t timeVariance[fIslNum];
519
520 //reset the "sets" functions
521 if (fIslNum <1)
522 fIsl->SetIslNum(0);
523
524 for(Int_t i = 0; i<20 ;i++)
525 {
526 fIsl->SetPixNum(i,-1);
527 fIsl->SetSigToNoise(i,-1);
528 fIsl->SetTimeSpread(i,-1);
529
530 for(Int_t idx = 0; idx<npix; idx++)
531 {
532 fIsl->SetIslId(idx, -1);
533 fIsl->SetArrivalTime(i, idx, -1 );
534 }
535 }
536
537
538 for(Int_t i = 1; i<=fIslNum ; i++)
539 {
540 Int_t n = 0;
541 Int_t ncore = 0;
542
543 Float_t MIN = 10000;
544 Float_t MAX = 0;
545
546 signal = 0;
547 noise = 0;
548 fPixNum[i-1] = 0;
549 timeVariance[i-1] = 0;
550
551 for(Int_t idx=0 ; idx<npix ; idx++)
552 {
553
554 MCerPhotPix *pix = fEvt->GetPixById(idx);
555 const MPedestalPix &ped = (*fPed)[idx];
556 const MArrivalTimePix &timepix = (*fTime)[idx];
557
558 if (pix == NULL) break;
559
560 if (vect[i][idx] ==1 && fEvt->IsPixelUsed(idx)){
561
562 fPixNum[i-1]++;
563 signal += pix->GetNumPhotons() * (fCam->GetPixRatio(idx));
564 noise += pow(ped.GetPedestalRms(),2);
565
566 time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
567
568
569 if (fEvt->IsPixelCore(idx)){
570
571 if (time[n] > MAX)
572 MAX = time[n];
573 if (time[n] < MIN)
574 MIN = time[n];
575
576 ncore++;
577 }
578
579 fIsl->SetIslId(idx, i-1);
580 fIsl->SetArrivalTime(i-1, idx, time[n]);
581
582 n++;
583 }
584
585 }
586
587 // Float_t mean = timeVariance[i-1]/ncore;
588
589 timeVariance[i-1] = 0;
590
591 timeVariance[i-1] = (MAX - MIN)/ncore;
592 timeVariance[i-1] = (MAX - MIN)/ncore;
593
594 fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
595
596 fIsl->SetPixNum(i-1,fPixNum[i-1]);
597 fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
598 fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
599
600 }
601
602 fIsl->SetReadyToSave();
603
604 return 1;
605}
Note: See TracBrowser for help on using the repository browser.