Changeset 5466
- Timestamp:
- 11/24/04 16:33:39 (20 years ago)
- Location:
- trunk/MagicSoft/Mars/mtemp/mifae
- Files:
-
- 3 edited
-
Changelog (modified) (1 diff)
-
library/MImgIsland.cc (modified) (2 diffs)
-
library/MIslandsCalc.cc (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mtemp/mifae/Changelog
r5379 r5466 18 18 19 19 -*-*- END OF LINE -*-*- 20 21 2004/11/24 Ester Aliu 22 * library/MIslandsCalc.h 23 - Remove a memmory leak, and add a sort function to order the 24 pixels by index after their interpolation 25 - Other minor changes 26 * library/MImgIsland.[h,cc] 27 - Reset some variables 20 28 21 29 2004/11/11 Ester Aliu -
trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc
r5379 r5466 37 37 fName = name ? name : "MImgIsland"; 38 38 fTitle = title ? title : "Storage container for one island information"; 39 39 // cout << "Creo 1 isla" << endl; 40 40 Reset(); 41 } 42 MImgIsland::~MImgIsland() 43 { 44 // cout << "Destruyo 1 isla" << endl; 41 45 } 42 46 … … 44 48 { 45 49 fPixNum = 0; 46 fSigToNoise = 0; 47 fTimeSpread = 0; 48 fMeanX = 0; 49 fMeanY = 0; 50 fDist = 0; 51 fDistW = 0; 52 fDistL = 0; 50 fSigToNoise = -1; 51 fTimeSpread = -1; 52 fMeanX = -1000; 53 fMeanY = -1000; 54 fDist = -1; 55 fDistW = -1; 56 fDistL = -1; 57 fDistS = -1; 58 59 fWidth = -1; 60 fLength = -1; 61 fSizeIsl = -1; 62 fAlpha = -1000; 53 63 54 64 fPixList.Reset(-1); 55 65 fPeakPulse.Reset(-1); 66 56 67 } 57 68 -
trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc
r5379 r5466 35 35 // - fPixNum // number of pixels in the island 36 36 // - fSigToNoise // signal to noise of the island 37 // - fTimeSpread // mean arrival time spreadof the island37 // - fTimeSpread // "time" of the island 38 38 // - fMeanX // mean X position of the island 39 39 // - fMeanY // mean Y position of the island 40 40 // - fDist // dist between an island and the continent 41 // - fLength // major axis of the largerisland ellipse42 // - fWidth // minor axis of the largerisland ellipse41 // - fLength // major axis of the island ellipse 42 // - fWidth // minor axis of the island ellipse 43 43 // - fDistL // dist divided by lenght of the larger island 44 44 // - fDistW // dist divided by width of the larger island … … 192 192 Float_t meanX[numisl]; 193 193 Float_t meanY[numisl]; 194 Float_t length[numisl]; 195 Float_t width[numisl]; 194 196 Float_t dist[numisl]; 195 197 Float_t distL[numisl]; … … 197 199 Float_t distS[numisl]; 198 200 199 Float_t size[numisl], sizeLargeIsl, length, width, distance, alpha, alphaW, sizetot; 200 201 Float_t X = 0; 202 Float_t Y = 0; 201 202 Float_t size[numisl], sizeLargeIsl, distance, alpha, alphaW, sizetot; 203 203 204 sizeLargeIsl = 0; 204 205 alphaW = 0; … … 211 212 212 213 imgIsl->InitSize(num[i]); 213 214 214 215 Int_t n = 0; 215 216 Float_t MIN = 10000.; 217 Float_t MAX = 0.; 218 216 217 Float_t minTime = 10000.; 218 Float_t maxTime = 0.; 219 220 Float_t minX = 10000.; 221 Float_t maxX = 0.; 222 223 Float_t minY = 10000.; 224 Float_t maxY = 0.; 225 219 226 signal = 0; 220 227 noise = 0; 221 228 222 229 size[i-1] = 0; 223 230 meanX[i-1] = 0; 224 231 meanY[i-1] = 0; 225 232 dist[i-1] = 0; 226 233 227 234 PixelNumIsl[i-1] = 0; 228 235 timeVariance[i-1] = 0; 229 230 236 231 237 for(Int_t idx=0 ; idx<nPix ; idx++) 232 238 { … … 247 253 if (i == 1) 248 254 sizeLargeIsl += nphot; 249 255 250 256 meanX[i-1] += nphot * gpix2.GetX(); 251 257 meanY[i-1] += nphot * gpix2.GetY(); 252 258 253 259 time[i-1] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain(); 254 260 255 261 imgIsl->SetPixList(PixelNumIsl[i-1]-1,pix->GetPixId()); 256 262 imgIsl->SetPeakPulse(PixelNumIsl[i-1]-1,time[i-1]); … … 259 265 if (fEvt->IsPixelCore(idx)){ 260 266 261 if (time[i-1] > MAX) 262 MAX = time[i-1]; 263 if (time[i-1] < MIN) 264 MIN = time[i-1]; 265 266 } 267 267 if (time[i-1] > maxTime){ 268 maxTime = time[i-1]; 269 maxX = gpix2.GetX(); 270 maxY = gpix2.GetY(); 271 } 272 if (time[i-1] < minTime){ 273 minTime = time[i-1]; 274 minX = gpix2.GetX(); 275 minY = gpix2.GetY(); 276 } 277 278 } 268 279 n++; 269 280 } 270 281 } 271 282 272 283 meanX[i-1] /= size[i-1]; 273 284 meanY[i-1] /= size[i-1]; 274 275 276 if (i == 1){ 277 X = meanX[i-1]; 278 Y = meanY[i-1]; 279 } 280 281 dist[i-1] = TMath::Power(meanX[i-1]-X,2) + TMath::Power(meanY[i-1]-Y,2); 285 286 dist[i-1] = TMath::Power(meanX[i-1]-meanX[0],2) + TMath::Power(meanY[i-1]-meanY[i-1],2); 282 287 dist[i-1] = TMath::Sqrt(dist[i-1]); 283 284 timeVariance[i-1] = MAX-MIN; 288 289 //timeVariance[i-1] = (maxTime-minTime); 290 291 if (maxX!=minX && maxY!=minY) 292 timeVariance[i-1] = (maxTime-minTime)/sqrt(TMath::Power(maxX-minX,2) + TMath::Power(maxY-minY,2)); 293 else 294 timeVariance[i-1] = -1; 285 295 286 296 //noise = 0, in the case of MC w/o noise 287 297 if (noise == 0) noise = 1; 288 298 289 299 SigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise); 290 300 291 301 imgIsl->SetPixNum(PixelNumIsl[i-1]); 292 302 imgIsl->SetSigToNoise(SigToNoise[i-1]); … … 296 306 imgIsl->SetDist(dist[i-1]); 297 307 imgIsl->SetSizeIsl(size[i-1]); 308 309 310 // sanity check: if one island has 2 or less pixels 311 if (num[i]>2){ 312 313 314 // calculate width and lenght of each island 315 316 Double_t corrxx=0; // [m^2] 317 Double_t corrxy=0; // [m^2] 318 Double_t corryy=0; // [m^2] 319 320 for(Int_t idx=0 ; idx<nPix ; idx++){ 321 322 MCerPhotPix *pix = fEvt->GetPixById(idx); 323 if(!pix) continue; 324 const MGeomPix &gpix3 = (*fCam)[pix->GetPixId()]; 325 const Float_t nphot = pix->GetNumPhotons(); 326 327 if (vect[i][idx]==1){ 328 329 const Float_t dx = gpix3.GetX() - meanX[i-1]; // [mm] 330 const Float_t dy = gpix3.GetY() - meanY[i-1]; // [mm] 331 332 corrxx += nphot * dx*dx; // [mm^2] 333 corrxy += nphot * dx*dy; // [mm^2] 334 corryy += nphot * dy*dy; // [mm^2] 335 336 } 337 } 338 339 const Double_t d0 = corryy - corrxx; 340 const Double_t d1 = corrxy*2; 341 const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1); 342 const Double_t tand = d2 / d1; 343 const Double_t tand2 = tand*tand; 344 345 const Double_t s2 = tand2+1; 346 const Double_t s = TMath::Sqrt(s2); 347 348 const Double_t CosDelta = 1.0/s; // need these in derived classes 349 const Double_t SinDelta = tand/s; // like MHillasExt 350 351 const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/size[i-1]; 352 const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/size[i-1]; 353 354 // 355 // fLength^2 is the second moment along the major axis of the ellipse 356 // fWidth^2 is the second moment along the minor axis of the ellipse 357 // 358 // From the algorithm we get: fWidth <= fLength is always true 359 // 360 // very small numbers can get negative by rounding 361 // 362 length[i-1] = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm] 363 width[i-1] = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm] 364 365 366 // alpha calculation 367 368 const Double_t mx = meanX[i-1]; // [mm] 369 const Double_t my = meanY[i-1]; // [mm] 370 371 //FIXME: xpos, ypos from MSrcPos 372 const Double_t xpos = 0.; 373 const Double_t ypos = 0.; 374 375 const Double_t sx = mx - xpos; // [mm] 376 const Double_t sy = my - ypos; // [mm] 377 378 const Double_t sd = SinDelta; // [1] 379 const Double_t cd = CosDelta; // [1] 380 381 // 382 // Distance from source position to center of ellipse. 383 // If the distance is 0 distance, Alpha is not specified. 384 // The calculation has failed and returnes kFALSE. 385 // 386 distance = TMath::Sqrt(sx*sx + sy*sy); // [mm] 387 if (distance==0){ 388 389 for (Int_t l = 0; l< nVect; l++) 390 delete [] vect[l]; 391 392 delete [] vect; 393 delete [] num; 394 395 return 1; 396 } 397 398 // 399 // Calculate Alpha and Cosda = cos(d,a) 400 // The sign of Cosda will be used for quantities containing 401 // a head-tail information 402 // 403 // *OLD* const Double_t arg = (sy-tand*sx) / (dist*sqrt(tand*tand+1)); 404 // *OLD* fAlpha = asin(arg)*kRad2Deg; 405 // 406 const Double_t arg2 = cd*sx + sd*sy; // [mm] 407 if (arg2==0){ 408 409 for (Int_t l = 0; l< nVect; l++) 410 delete [] vect[l]; 411 412 delete [] vect; 413 delete [] num; 414 415 return 2; 416 } 417 418 const Double_t arg1 = cd*sy - sd*sx; // [mm] 419 420 // 421 // Due to numerical uncertanties in the calculation of the 422 // square root (dist) and arg1 it can happen (in less than 1e-5 cases) 423 // that the absolute value of arg exceeds 1. Because this uncertainty 424 // results in an Delta Alpha which is still less than 1e-3 we don't care 425 // about this uncertainty in general and simply set values which exceed 426 // to 1 saving its sign. 427 // 428 const Double_t arg = arg1/distance; 429 alpha = TMath::Abs(arg)>1 ? TMath::Sign(90., arg) : TMath::ASin(arg)*TMath::RadToDeg(); // [deg] 430 431 alphaW += alpha*size[i-1]; 432 sizetot += size[i-1]; 433 434 //////////////////////////////////////// 435 436 distL[i-1]=dist[i-1]/length[0]; 437 distW[i-1]=dist[i-1]/width[0]; 438 distS[i-1]= dist[i-1]/size[0]; 439 440 imgIsl->SetLength(length[i-1]); 441 imgIsl->SetWidth(width[i-1]); 442 443 imgIsl->SetDistL(distL[i-1]); 444 imgIsl->SetDistW(distW[i-1]); 445 imgIsl->SetDistS(distS[i-1]); 446 447 imgIsl->SetAlpha(alpha); 448 449 } 298 450 299 451 fIsl->GetList()->Add(imgIsl); 300 452 301 453 } 302 303 304 //Length and Width of the larger island according the definition of the hillas parameters 305 306 // calculate 2nd moments 307 // --------------------- 308 Double_t corrxx=0; // [m^2] 309 Double_t corrxy=0; // [m^2] 310 Double_t corryy=0; // [m^2] 311 312 for(Int_t idx=0 ; idx<nPix ; idx++) 313 { 314 MCerPhotPix *pix = fEvt->GetPixById(idx); 315 if(!pix) continue; 316 const MGeomPix &gpix3 = (*fCam)[pix->GetPixId()]; 317 const Float_t nphot = pix->GetNumPhotons(); 318 319 // if (pix == NULL) break; 320 321 if (vect[1][idx]==1){ 322 323 const Float_t dx = gpix3.GetX() - X; // [mm] 324 const Float_t dy = gpix3.GetY() - Y; // [mm] 325 326 327 corrxx += nphot * dx*dx; // [mm^2] 328 corrxy += nphot * dx*dy; // [mm^2] 329 corryy += nphot * dy*dy; // [mm^2] 330 331 } 332 } 333 334 335 // calculate the hillas parameters Width and Length 336 337 MImgIsland *imgIsl = new MImgIsland; 338 TIter Next(fIsl->GetList()); 339 //Next.Reset(); 340 341 Int_t i = 1; 342 while ((imgIsl=(MImgIsland*)Next())) { 343 344 const Double_t d0 = corryy - corrxx; 345 const Double_t d1 = corrxy*2; 346 const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1); 347 const Double_t tand = d2 / d1; 348 const Double_t tand2 = tand*tand; 349 350 const Double_t s2 = tand2+1; 351 const Double_t s = TMath::Sqrt(s2); 352 353 const Double_t CosDelta = 1.0/s; // need these in derived classes 354 const Double_t SinDelta = tand/s; // like MHillasExt 355 356 const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/size[i-1]; 357 const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/size[i-1]; 358 359 // 360 // fLength^2 is the second moment along the major axis of the ellipse 361 // fWidth^2 is the second moment along the minor axis of the ellipse 362 // 363 // From the algorithm we get: fWidth <= fLength is always true 364 // 365 // very small numbers can get negative by rounding 366 // 367 length = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm] 368 width = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm] 369 370 371 // alpha calculation 372 373 const Double_t mx = imgIsl->GetMeanX(); // [mm] 374 const Double_t my = imgIsl->GetMeanY(); // [mm] 375 376 //FIXME: xpos, ypos from MSrcPos 377 const Double_t xpos = 0.; 378 const Double_t ypos = 0.; 379 380 const Double_t sx = mx - xpos; // [mm] 381 const Double_t sy = my - ypos; // [mm] 382 383 const Double_t sd = SinDelta; // [1] 384 const Double_t cd = CosDelta; // [1] 385 386 // 387 // Distance from source position to center of ellipse. 388 // If the distance is 0 distance, Alpha is not specified. 389 // The calculation has failed and returnes kFALSE. 390 // 391 distance = TMath::Sqrt(sx*sx + sy*sy); // [mm] 392 if (distance==0) 393 return 1; 394 395 // 396 // Calculate Alpha and Cosda = cos(d,a) 397 // The sign of Cosda will be used for quantities containing 398 // a head-tail information 399 // 400 // *OLD* const Double_t arg = (sy-tand*sx) / (dist*sqrt(tand*tand+1)); 401 // *OLD* fAlpha = asin(arg)*kRad2Deg; 402 // 403 const Double_t arg2 = cd*sx + sd*sy; // [mm] 404 if (arg2==0) 405 return 2; 406 407 const Double_t arg1 = cd*sy - sd*sx; // [mm] 408 409 // 410 // Due to numerical uncertanties in the calculation of the 411 // square root (dist) and arg1 it can happen (in less than 1e-5 cases) 412 // that the absolute value of arg exceeds 1. Because this uncertainty 413 // results in an Delta Alpha which is still less than 1e-3 we don't care 414 // about this uncertainty in general and simply set values which exceed 415 // to 1 saving its sign. 416 // 417 const Double_t arg = arg1/distance; 418 alpha = TMath::Abs(arg)>1 ? TMath::Sign(90., arg) : TMath::ASin(arg)*TMath::RadToDeg(); // [deg] 419 420 alphaW += alpha*size[i-1]; 421 sizetot += size[i-1]; 422 423 //////////////////////////////////////// 424 425 distL[i-1]=dist[i-1]/length; 426 distW[i-1]=dist[i-1]/width; 427 distS[i-1]= dist[i-1]/size[i-1]; 428 429 imgIsl->SetLength(length); 430 imgIsl->SetWidth(width); 431 432 imgIsl->SetDistL(distL[i-1]); 433 imgIsl->SetDistW(distW[i-1]); 434 imgIsl->SetDistS(distS[i-1]); 435 436 imgIsl->SetAlpha(alpha); 437 i++; 438 } 439 440 454 441 455 fIsl->SetAlphaW(alphaW/sizetot); 442 456 //fIsl->SetReadyToSave(); 443 444 for (Int_t i = 0; i< nVect; i++) 445 delete [] vect[i]; 446 447 delete vect; 457 458 for (Int_t l = 0; l< nVect; l++) 459 delete [] vect[l]; 460 461 delete [] vect; 462 delete [] num; 448 463 449 464 return kTRUE; … … 481 496 MCerPhotPix *pix; 482 497 483 //loop over all pixels 484 MCerPhotEvtIter Next(fEvt, kFALSE); 485 498 // Loop over used pixels only 499 TIter Next(*fEvt); 500 501 //after the interpolation of the pixels, these can be disordered by index. This is important for this algorithm of calculating islands 502 fEvt->Sort(); 503 486 504 while ((pix=static_cast<MCerPhotPix*>(Next()))) 487 505 { … … 493 511 if( fEvt->IsPixelUsed(idx)) 494 512 { 495 //cout << 513 //cout <<idx <<endl; 496 514 sflag = 0; 497 515 … … 650 668 MCerPhotPix *pix; 651 669 652 //1st loop over all pixels 653 MCerPhotEvtIter Next0(fEvt, kFALSE); 670 //after the interpolation of the pixels, these can be disordered by index. This is important for this algorithm of calculating islands 671 fEvt->Sort(); 672 673 // 1st loop over used pixels only 674 TIter Next0(*fEvt); 654 675 655 676 while ((pix=static_cast<MCerPhotPix*>(Next0()))) … … 671 692 } 672 693 694 673 695 //2nd loop over all pixels 674 MCerPhotEvtIter Next(fEvt, kFALSE);675 696 TIter Next(*fEvt); 697 676 698 while ((pix=static_cast<MCerPhotPix*>(Next()))) 677 699 { … … 713 735 numisl = nvect; 714 736 715 716 // Repeated Chain Corrections 717 718 Int_t ii, jj; 719 720 for(Int_t i = 1; i <= nvect; i++) { 721 for(Int_t j = 1; j <= nvect; j++) { 722 if (i!=j && zeros[j]!=1){ 723 control = 0; 724 for(Int_t k = 0; k < npix; k++) { 725 if (vect[i][k] == 1 && vect[j][k] == 1) { 726 control = 1; 727 break; 728 } 729 } 730 if(i<j) { 731 ii=i; 732 jj=j; 733 } 734 else{ 735 ii=j; 736 jj=i; 737 } 738 if (control == 1) { 739 for(Int_t k = 0; k < npix; k++) { 740 if(vect[jj][k] == 1) vect[ii][k] = 1; 741 vect[jj][k] = 0; 742 zeros[jj] = 1; 743 } 744 numisl = numisl-1; 745 } 746 } 747 } 748 } 737 // Repeated Chain Corrections 738 739 Int_t jmin = 0; 740 741 for(Int_t i = 1; i <= nvect; i++){ 742 control=0; 743 for(Int_t j = i+1; j <= nvect; j++){ 744 control = 0; 745 for(Int_t k = 0; k < npix; k++){ 746 if (vect[i][k] == 1 && vect[j][k] == 1){ 747 control = 1; 748 k=npix; 749 } 750 } 751 if (control == 1){ 752 for(Int_t k = 0; k < npix; k++){ 753 if(vect[j][k] == 1) vect[i][k] = 1; 754 vect[j][k] = 0; 755 zeros[j] = 1; 756 } 757 numisl = numisl-1; 758 } 759 } 760 761 for(Int_t j = 1; j <= i-1; j++){ 762 for(Int_t k = 0; k < npix; k++){ 763 if (vect[i][k] == 1 && vect[j][k] == 1){ 764 control = 2; 765 jmin=j; 766 k=npix; 767 j=i; 768 } 769 } 770 771 if (control == 2){ 772 for (Int_t k = 0; k < npix; k++){ 773 if(vect[i][k]==1) vect[jmin][k]=1; 774 vect[i][k] = 0; 775 zeros[i] = 1; 776 } 777 numisl = numisl-1; 778 } 779 } 780 } 749 781 750 782 Int_t l = 1;
Note:
See TracChangeset
for help on using the changeset viewer.
