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

Last change on this file since 3976 was 3976, checked in by aliu, 21 years ago
*** empty log message ***
File size: 8.9 KB
Line 
1/* ======================================================================== *\
2!
3!
4! Author(s): Ester Aliu, 2/2004 <aliu@ifae.es>
5!
6!
7\* ======================================================================== */
8
9/////////////////////////////////////////////////////////////////////////////
10//
11// MIslandCalc
12//
13// The Island Calc task calculates the some islands parameters for each
14// of the events such as:
15//
16// - fIslNum // number of islands
17// - fIslId[577] // island Id
18// - fPixNum[20]; // number of pixels in the island
19// - fSigToNoise[20]; // signal to noise of the island
20// - fTimeHi[20][577] // mean of the arrival time in the hi gain
21// - fTimeLo[20][577] // mean of the arrival time in the lo gain
22// - fTimeSpreadHi[20]; // mean arrival time spread of the island
23// - fTimeSpreadLo[20]; // mean arrival time spread of the island
24//
25//
26// Input Containers:
27// MGeomCam
28// MCerPhotEvt
29// MPedestalCam
30// MArrivalTimeCam
31//
32// Output Containers:
33// MIslands
34//
35/////////////////////////////////////////////////////////////////////////////
36#include "MIslandCalc.h"
37
38#include <stdlib.h> // atof
39#include <fstream> // ofstream, SavePrimitive
40
41#include "MLog.h"
42#include "MLogManip.h"
43
44#include "MIslands.h"
45
46#include "MParList.h"
47
48#include "MGeomPix.h"
49#include "MGeomCam.h"
50
51#include "MCerPhotPix.h"
52#include "MCerPhotEvt.h"
53
54#include "MPedestalCam.h"
55#include "MPedestalPix.h"
56
57#include "MArrivalTimeCam.h"
58#include "MArrivalTimePix.h"
59
60ClassImp(MIslandCalc);
61
62
63using namespace std;
64
65// --------------------------------------------------------------------------
66//
67// Default constructor.
68//
69MIslandCalc::MIslandCalc(const char* name, const char* title)
70 : fIsl(NULL)
71{
72 fName = name ? name : "MIslandCalc";
73 fTitle = title ? title : "Calculate island parameters";
74}
75
76
77// --------------------------------------------------------------------------
78Int_t MIslandCalc::PreProcess (MParList *pList)
79{
80 fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
81 if (!fCam)
82 {
83 *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
84 return kFALSE;
85 }
86
87 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
88 if (!fEvt)
89 {
90 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
91 return kFALSE;
92 }
93
94 fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
95 if (!fPed)
96 {
97 *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
98 return kFALSE;
99 }
100
101 fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam"));
102 if (!fTime)
103 {
104 *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl;
105 return kFALSE;
106 }
107
108 if (strlen(fIslName) > 0)
109 {
110 fIsl = (MIslands*)pList->FindCreateObj("MIslands", AddSerialNumber(fIslName));
111 cout << "kk1" << endl;
112 }
113 else
114 {
115 fIsl = (MIslands*)pList->FindCreateObj(AddSerialNumber("MIslands"));
116 cout << "kk2" << endl;
117 }
118 if (!fIsl)
119 return kFALSE;
120
121 return kTRUE;
122}
123
124
125Int_t MIslandCalc::Process()
126{
127
128 Float_t noise;
129 Float_t signal;
130
131 Int_t npix = 577;
132 Int_t sflag;
133 Int_t control;
134
135 Int_t nvect = 0;
136 Int_t fIslNum = 0;
137
138 Int_t vect[50][577];
139
140 Int_t zeros[50];
141
142 for(Int_t m = 0; m < 50 ; m++)
143 for(Int_t n = 0; n < npix ; n++)
144 vect[m][n] = 0;
145
146 for(Int_t n = 0; n < 50 ; n++)
147 zeros[n] = 0;
148
149
150 //loop over all pixels
151 for(Int_t idx=0 ; idx<npix ; idx++)
152 {
153 const MGeomPix &gpix = (*fCam)[idx];
154 const Int_t nnmax = gpix.GetNumNeighbors();
155
156 if(idx<0 || idx==npix)
157 {
158 cout << "Pixel (software) index out of range [0-576]. Available number of pixels=" << npix << endl;
159 return 1;
160 }
161
162
163 if( fEvt->IsPixelUsed(idx))
164 {
165 sflag = 0;
166
167 for(Int_t j=0; j < nnmax ; j++)
168 {
169 const Int_t idx2 = gpix.GetNeighbor(j);
170
171 if (idx2 < idx)
172 {
173 for(Int_t k = 1; k <= nvect; k++)
174 {
175 if (vect[k][idx2] == 1)
176 {
177 sflag = 1;
178 vect[k][idx] = 1;
179 }
180 }
181 }
182 }
183
184 if (sflag == 0)
185 {
186 nvect++;
187 vect[nvect][idx] = 1;
188 }
189
190 }
191 }
192
193 fIslNum = nvect;
194
195
196 // Repeated Chain Corrections
197
198 for(Int_t i = 1; i <= nvect; i++)
199 {
200 for(Int_t j = i+1; j <= nvect; j++)
201 {
202 control = 0;
203 for(Int_t k = 0; k < npix; k++)
204 {
205
206 if (vect[i][k] == 1 && vect[j][k] == 1)
207 {
208 control = 1;
209 break;
210 }
211 }
212 if (control == 1)
213 {
214 for(Int_t k = 0; k < npix; k++)
215 {
216 if(vect[j][k] == 1)
217 vect[i][k] = 1;
218 vect[j][k] = 0;
219 zeros[j] = 1;
220 }
221 fIslNum = fIslNum-1;
222 }
223
224 }
225 }
226
227
228 Int_t l = 1;
229
230 for(Int_t i = 1; i<= nvect ; i++)
231 {
232 if (zeros[i] == 0)
233 {
234 for(Int_t k = 0; k<npix; k++)
235 {
236 vect[l][k] = vect[i][k];
237 }
238 l++;
239 }
240 }
241
242
243 //set the number of islands in one event
244 fIsl->SetIslNum(fIslNum);
245
246 //examine each island...
247 Int_t fPixNum[fIslNum];
248 Float_t fSigToNoise[fIslNum];
249 Float_t time[577];
250 Float_t timeVariance[fIslNum];
251 //Float_t timeHi[577], timeLo[577];
252 //Float_t timeVarianceHi[fIslNum];
253 //Float_t timeVarianceLo[fIslNum];
254
255 //reset the "sets" functions
256 if (fIslNum <1)
257 fIsl->SetIslNum(0);
258
259 for(Int_t i = 0; i<20 ;i++)
260 {
261 fIsl->SetPixNum(i,-1);
262 fIsl->SetSigToNoise(i,-1);
263 fIsl->SetTimeSpread(i,-1);
264 // fIsl->SetTimeSpreadHi(i,-1);
265 // fIsl->SetTimeSpreadLo(i,-1);
266
267 for(Int_t idx = 0; idx<npix; idx++)
268 {
269 fIsl->SetIslId(idx, -1);
270 fIsl->SetArrivalTime(i, idx, -1 );
271 // fIsl->SetArrivalTimeHiGain(i, idx, -1 );
272 // fIsl->SetArrivalTimeLoGain(i, idx, -1);
273 }
274 }
275
276
277 for(Int_t i = 1; i<=fIslNum ; i++)
278 {
279 Int_t n = 0;
280 Int_t ncore = 0;
281
282 Float_t MINhi = 10000;
283 Float_t MINlo = 10000;
284 Float_t MAXhi = 0;
285 Float_t MAXlo = 0;
286 Float_t MIN = 10000;
287 Float_t MAX = 0;
288
289 //cout << "Isl #" << i << endl;
290
291 signal = 0;
292 noise = 0;
293 fPixNum[i-1] = 0;
294 timeVariance[i-1] = 0;
295 // timeVarianceHi[i-1] = 0;
296 // timeVarianceLo[i-1] = 0;
297
298 for(Int_t idx=0 ; idx<npix ; idx++)
299 {
300 const MCerPhotPix &pix = (*fEvt)[idx];
301 const MPedestalPix &ped = (*fPed)[idx];
302 const MArrivalTimePix &timepix = (*fTime)[idx];
303
304 if (vect[i][idx]==1){
305
306 fPixNum[i-1]++;
307 signal += pix.GetNumPhotons() * (fCam->GetPixRatio(idx));
308 noise += pow(ped.GetPedestalRms(),2);
309
310 time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
311
312 // timeHi[n] = timepix.GetArrivalTimeHiGain();
313 // cout << "timeHi" << timeHi[n] << endl;
314
315 //timeLo[n] = timepix.GetArrivalTimeLoGain();
316 //cout << "timeHi" << timeLo[n] << endl;
317
318 if (fEvt->IsPixelCore(idx)){
319
320 /* if (timeHi[n] < MINhi)
321 MINhi = timeHi[n];
322 if (timeLo[n] < MINlo)
323 MINlo = timeLo[n];
324
325 if (timeHi[n] > MAXhi)
326 MAXhi = timeHi[n];
327 if (timeLo[n] > MAXlo)
328 MAXlo = timeLo[n];
329
330 */
331
332 if (time[n] > MAX)
333 MAX = time[n];
334 if (time[n] < MIN)
335 MIN = time[n];
336
337 //control2[n] = 1;
338 //timeVarianceHi[i-1] += timeHi[n];
339 //timeVarianceLo[i-1] += timeLo[n];
340 ncore++;
341 }
342
343 fIsl->SetIslId(idx, i-1);
344 fIsl->SetArrivalTime(i-1, idx, time[n]);
345 // fIsl->SetArrivalTimeHiGain(i-1, idx, timeHi[n]);
346 // fIsl->SetArrivalTimeLoGain(i-1, idx, timeLo[n]);
347
348 n++;
349 }
350
351 }
352
353 Float_t mean = timeVariance[i-1]/ncore;
354 // Float_t meanHi = timeVarianceHi[i-1]/ncore;
355 // Float_t meanLo = timeVarianceLo[i-1]/ncore;
356
357 timeVariance[i-1] = 0;
358 //timeVarianceHi[i-1] = 0;
359 //timeVarianceLo[i-1] = 0;
360
361 /*
362 //old
363 for (Int_t k = 0; k <n ; k++)
364 {
365 if (control2[k] == 1){
366
367 timeVarianceHi[i-1] += pow(timeHi[k]- meanHi,2);
368 timeVarianceLo[i-1] += pow(timeLo[k]- meanLo,2);
369 }
370 }
371 timeVarianceHi[i-1] = sqrt(timeVarianceHi[i-1]/(ncore-1));
372 timeVarianceLo[i-1] = sqrt(timeVarianceLo[i-1]/(ncore-1)); */
373
374 //timeVarianceHi[i-1] = (MAXhi - MINhi)/ncore;
375 //timeVarianceLo[i-1] = (MAXlo - MINlo)/ncore;
376
377 timeVariance[i-1] = (MAX - MIN)/ncore;
378 timeVariance[i-1] = (MAX - MIN)/ncore;
379
380 fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
381
382 fIsl->SetPixNum(i-1,fPixNum[i-1]);
383 fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
384 fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
385 //fIsl->SetTimeSpreadHi(i-1, timeVarianceHi[i-1]);
386 //fIsl->SetTimeSpreadLo(i-1, timeVarianceLo[i-1]);
387
388 //cout << " TimeHi Spread: " << timeVarianceHi[i-1] << endl;
389 //cout << " # core pixels: " << ncore << endl;
390 //cout << " # used pixels: " << n << endl;
391 //cout << " SigToNoise: " << fSigToNoise[i-1] << endl << endl;
392
393 }
394
395 fIsl->SetReadyToSave();
396
397 return 1;
398
399}
400
401
Note: See TracBrowser for help on using the repository browser.