| 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 |
|
|---|
| 60 | ClassImp(MIslandCalc);
|
|---|
| 61 |
|
|---|
| 62 |
|
|---|
| 63 | using namespace std;
|
|---|
| 64 |
|
|---|
| 65 | // --------------------------------------------------------------------------
|
|---|
| 66 | //
|
|---|
| 67 | // Default constructor.
|
|---|
| 68 | //
|
|---|
| 69 | MIslandCalc::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 | // --------------------------------------------------------------------------
|
|---|
| 78 | Int_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 |
|
|---|
| 125 | Int_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 |
|
|---|
