Changeset 5379
- Timestamp:
- 11/11/04 09:55:25 (20 years ago)
- Location:
- trunk/MagicSoft/Mars/mtemp/mifae
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mtemp/mifae/Changelog
r5295 r5379 18 18 19 19 -*-*- END OF LINE -*-*- 20 21 2004/11/11 Ester Aliu 22 * library/MIslands.h 23 - Add the variable AlphaW 24 * library/MImgIsland.[h,cc] 25 - Add alpha for each island 26 * library/MIslandsCalc.cc 27 - Bug solved 28 * library/MIslandsClean.cc 29 - More different island cleanings added 20 30 21 31 2004/10/20 Javier Rico -
trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc
r5170 r5379 16 16 ! 17 17 ! Author(s): Javier Rico 04/2004 <mailto:jrico@ifae.es> 18 ! 18 ! Ester Aliu 10/2004 <mailto:aliu@ifae.es> (update according 19 ! the new classes of the islands) 19 20 ! Copyright: MAGIC Software Development, 2000-2004 20 21 ! -
trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc
r5170 r5379 54 54 fPixList.Reset(-1); 55 55 fPeakPulse.Reset(-1); 56 57 56 } 58 57 … … 105 104 *fLog << inf << " DistW = " << fDistW << endl; 106 105 *fLog << inf << " DistS = " << fDistS << endl; 106 *fLog << inf << " Alpha = " << fAlpha << endl; 107 107 108 108 } -
trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h
r5170 r5379 22 22 Float_t fSigToNoise; 23 23 Float_t fTimeSpread; 24 Float_t fMeanX;25 Float_t fMeanY;26 Float_t fWidth;27 Float_t fLength;28 24 Float_t fDist; 29 25 Float_t fDistL; 30 26 Float_t fDistW; 31 27 Float_t fDistS; 28 29 Float_t fWidth; 30 Float_t fLength; 31 Float_t fSizeIsl; 32 Float_t fMeanX; 33 Float_t fMeanY; 34 35 Float_t fAlpha; 36 32 37 33 38 TArrayI fPixList; … … 42 47 Float_t GetSigToNoise() { return fSigToNoise; } 43 48 Float_t GetTimeSpread() { return fTimeSpread; } 49 Float_t GetDist() { return fDist; } 50 Float_t GetDistL() { return fDistL; } 51 Float_t GetDistW() { return fDistW; } 52 Float_t GetDistS() { return fDistS; } 53 54 //hillas parameters 55 Float_t GetSizeIsl() { return fSizeIsl; } 44 56 Float_t GetMeanX() { return fMeanX; } 45 57 Float_t GetMeanY() { return fMeanY; } 46 58 Float_t GetWidth() { return fWidth; } 47 59 Float_t GetLength() { return fLength; } 48 Float_t GetDist() { return fDist; }49 Float_t GetDistL() { return fDistL; }50 Float_t GetDistW() { return fDistW; }51 Float_t GetDistS() { return fDistS; }52 60 61 // hillas src parameters 62 Float_t GetAlpha() { return fAlpha; } 63 53 64 void InitSize(Int_t i); 54 65 UInt_t GetSize() const { return fPixList.GetSize(); } … … 59 70 void Reset(); 60 71 61 void SetPixNum (Int_t i) { fPixNum = i;} 62 void SetSigToNoise(Float_t val) { fSigToNoise = val;} 63 void SetTimeSpread(Float_t val) { fTimeSpread = val;} 64 void SetMeanX (Float_t val) { fMeanX = val;} 65 void SetMeanY (Float_t val) { fMeanY = val;} 66 void SetDist (Float_t val) { fDist = val;} 67 void SetWidth (Float_t val) { fWidth = val;} 68 void SetLength (Float_t val) { fLength = val;} 69 void SetDistL (Float_t val) { fDistL = val;} 70 void SetDistW (Float_t val) { fDistW = val;} 71 void SetDistS (Float_t val) { fDistS = val;} 72 void SetPixNum (Int_t i) { fPixNum = i; } 73 void SetSigToNoise(Float_t val) { fSigToNoise = val; } 74 void SetTimeSpread(Float_t val) { fTimeSpread = val; } 75 void SetDist (Float_t val) { fDist = val; } 76 void SetDistL (Float_t val) { fDistL = val; } 77 void SetDistW (Float_t val) { fDistW = val; } 78 void SetDistS (Float_t val) { fDistS = val; } 79 80 //hillas parameters 81 void SetSizeIsl (Float_t val) { fSizeIsl = val; } 82 void SetMeanX (Float_t val) { fMeanX = val; } 83 void SetMeanY (Float_t val) { fMeanY = val; } 84 void SetWidth (Float_t val) { fWidth = val; } 85 void SetLength (Float_t val) { fLength = val; } 86 87 // hillas src parameters 88 void SetAlpha (Float_t val) { fAlpha = val; } 72 89 73 void SetPixList( const Int_t i, const Int_t id);74 void SetPeakPulse( const Int_t i, const Float_t time);90 void SetPixList( const Int_t i, const Int_t id); 91 void SetPeakPulse( const Int_t i, const Float_t time); 75 92 76 93 // void Paint(Option_t *opt=NULL); 77 94 void Print(Option_t *opt=NULL) const; 78 95 79 ClassDef(MImgIsland, 1) // Container that holds the island information96 ClassDef(MImgIsland, 2) // Container that holds the island information 80 97 81 98 }; -
trunk/MagicSoft/Mars/mtemp/mifae/library/MIslands.h
r5186 r5379 23 23 24 24 Int_t fIslNum; 25 Float_t fAlphaW; 25 26 26 27 public: … … 34 35 TList* GetList() const { return fIslands; } 35 36 UInt_t GetIslNum() const { return fIslands->GetSize(); } //number of islands 37 Float_t GetAlphaW() const { return fAlphaW; } //alpha weighted 36 38 37 void SetIslNum (Int_t i) { fIslNum = i; }38 39 void SetIslNum (Int_t i) { fIslNum = i; } 40 void SetAlphaW(Float_t val) { fAlphaW = val; } 39 41 // MHCamera& GetDisplay() { return fDisplay; } 40 42 -
trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc
r5170 r5379 180 180 fIsl->SetIslNum(numisl); 181 181 182 //cout << "code numisl " << numisl << endl; 182 183 //examine each island... 183 184 … … 196 197 Float_t distS[numisl]; 197 198 198 Float_t size[numisl], sizeLargeIsl, length, width ;199 Float_t size[numisl], sizeLargeIsl, length, width, distance, alpha, alphaW, sizetot; 199 200 200 201 Float_t X = 0; 201 202 Float_t Y = 0; 202 203 sizeLargeIsl = 0; 204 alphaW = 0; 205 sizetot = 0; 203 206 204 207 for(Int_t i = 1; i<=numisl ; i++) … … 292 295 imgIsl->SetMeanY(meanY[i-1]); 293 296 imgIsl->SetDist(dist[i-1]); 297 imgIsl->SetSizeIsl(size[i-1]); 294 298 295 299 fIsl->GetList()->Add(imgIsl); … … 333 337 MImgIsland *imgIsl = new MImgIsland; 334 338 TIter Next(fIsl->GetList()); 335 339 //Next.Reset(); 340 336 341 Int_t i = 1; 337 342 while ((imgIsl=(MImgIsland*)Next())) { … … 344 349 345 350 const Double_t s2 = tand2+1; 346 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 347 356 const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/size[i-1]; 348 357 const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/size[i-1]; … … 359 368 width = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm] 360 369 361 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 362 425 distL[i-1]=dist[i-1]/length; 363 426 distW[i-1]=dist[i-1]/width; 364 427 distS[i-1]= dist[i-1]/size[i-1]; 428 365 429 imgIsl->SetLength(length); 366 430 imgIsl->SetWidth(width); 431 367 432 imgIsl->SetDistL(distL[i-1]); 368 433 imgIsl->SetDistW(distW[i-1]); 369 434 imgIsl->SetDistS(distS[i-1]); 435 436 imgIsl->SetAlpha(alpha); 370 437 i++; 371 438 } 372 439 440 441 fIsl->SetAlphaW(alphaW/sizetot); 373 442 //fIsl->SetReadyToSave(); 374 443 … … 394 463 395 464 Int_t sflag; 396 Int_t control ;465 Int_t control = 0; 397 466 398 467 Int_t nvect = 0; … … 455 524 numisl = nvect; 456 525 457 526 458 527 // Repeated Chain Corrections 459 528 460 529 Int_t jmin = 0; 461 530 462 for(Int_t i = 1; i <= nvect; i++) 463 { 464 for(Int_t j = i+1; j <= nvect; j++) 465 { 466 control = 0; 467 for(Int_t k = 0; k < npix; k++) 468 { 469 if (vect[i][k] == 1 && vect[j][k] == 1) 470 { 471 control = 1; 472 break; 473 } 474 else if(vect[i][k] == 1 && vect[j][k] == 0){ 475 476 for(Int_t jj=1 ;jj<i ; jj++){ 477 478 if(vect[jj][k]==1){ 479 jmin = jj; 480 control = 2; 481 break; 482 } 483 } 484 } 485 } 486 487 488 if (control == 1) 489 { 490 for(Int_t k = 0; k < npix; k++) 491 { 492 if(vect[j][k] == 1) 493 vect[i][k] = 1; 494 vect[j][k] = 0; 495 zeros[j] = 1; 496 } 497 numisl = numisl-1; 498 } 499 500 if (control == 2) 501 { 502 for (Int_t k = 0; k < npix; k++) 503 { 504 if(vect[i][k]==1) 505 vect[jmin][k]=1; 506 vect[i][k] = 0; 507 zeros[i] = 1; 508 } 509 numisl = numisl-1; 510 } 511 } 512 } 531 for(Int_t i = 1; i <= nvect; i++){ 532 control=0; 533 for(Int_t j = i+1; j <= nvect; j++){ 534 control = 0; 535 for(Int_t k = 0; k < npix; k++){ 536 if (vect[i][k] == 1 && vect[j][k] == 1){ 537 control = 1; 538 k=npix; 539 } 540 } 541 if (control == 1){ 542 for(Int_t k = 0; k < npix; k++){ 543 if(vect[j][k] == 1) vect[i][k] = 1; 544 vect[j][k] = 0; 545 zeros[j] = 1; 546 } 547 numisl = numisl-1; 548 } 549 } 550 551 for(Int_t j = 1; j <= i-1; j++){ 552 for(Int_t k = 0; k < npix; k++){ 553 if (vect[i][k] == 1 && vect[j][k] == 1){ 554 control = 2; 555 jmin=j; 556 k=npix; 557 j=i; 558 } 559 } 560 561 if (control == 2){ 562 for (Int_t k = 0; k < npix; k++){ 563 if(vect[i][k]==1) vect[jmin][k]=1; 564 vect[i][k] = 0; 565 zeros[i] = 1; 566 } 567 numisl = numisl-1; 568 } 569 } 570 } 513 571 514 572 Int_t pixMAX = 0; … … 540 598 } 541 599 } 542 600 601 543 602 //the larger island will correspond to the 1st component of the vector 544 603 … … 657 716 // Repeated Chain Corrections 658 717 659 Int_t jmin = 0; 660 661 for(Int_t i = 1; i <= nvect; i++) 662 { 663 for(Int_t j = i+1; j <= nvect; j++) 664 { 665 control = 0; 666 for(Int_t k = 0; k < npix; k++) 667 { 668 669 if (vect[i][k] == 1 && vect[j][k] == 1) 670 { 671 control = 1; 672 break; 673 } 674 else if(vect[i][k] == 1 && vect[j][k] == 0){ 675 676 for(Int_t jj=1 ;jj<i ; jj++){ 677 678 if(vect[jj][k]==1){ 679 jmin = jj; 680 control = 2; 681 break; 682 } 683 } 684 } 685 } 686 687 688 if (control == 1) 689 { 690 for(Int_t k = 0; k < npix; k++) 691 { 692 if(vect[j][k] == 1) 693 vect[i][k] = 1; 694 vect[j][k] = 0; 695 zeros[j] = 1; 696 } 697 numisl = numisl-1; 698 } 699 700 if (control == 2) 701 { 702 for (Int_t k = 0; k < npix; k++) 703 { 704 if(vect[i][k]==1) 705 vect[jmin][k]=1; 706 vect[i][k] = 0; 707 zeros[i] = 1; 708 } 709 numisl = numisl-1; 710 } 711 712 } 713 } 714 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 } 715 749 716 750 Int_t l = 1; -
trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.cc
r5170 r5379 144 144 Int_t MIslandsClean::Process() 145 145 { 146 //147 //eliminates the island with a signal-to-noise148 //lower than a given limit149 //150 //if ( fIslandCleaningMethod == kNoTiming ){151 146 152 147 MImgIsland *imgIsl = new MImgIsland; … … 156 151 Int_t idPix = -1; 157 152 158 if ( fIslandCleaningMethod == 1 ) { 153 //////////////////////////////////////////////////// 154 // 155 // TIME SPREAD ISLAND CLEANING 156 // eliminates the island with a time spread 157 // higher than a given limit 158 /////////////////////////////////////////////////// 159 if( fIslandCleaningMethod == 0 ){ 160 while ((imgIsl=(MImgIsland*)Next())) { 161 162 //fIslCleanThreshold has different values, FIXME, put two variables 163 if(imgIsl->GetTimeSpread() > fIslCleanThres) 164 { 165 pixNum = imgIsl->GetPixNum(); 166 167 for(Int_t k = 0; k<pixNum; k++) 168 { 169 idPix = imgIsl->GetPixList(k); 170 MCerPhotPix &pix = (*fEvt)[idPix]; 171 pix.SetPixelUnused(); 172 } 173 } 174 } 175 176 } 177 178 179 //////////////////////////////////////////////////// 180 // 181 // SIGNAL TO NOISE ISLAND CLEANING 182 // eliminates the island with a signal-to-noise 183 // lower than a given limit 184 /////////////////////////////////////////////////// 185 else if ( fIslandCleaningMethod == 1 ) { 159 186 while ((imgIsl=(MImgIsland*)Next())) { 160 187 … … 173 200 } 174 201 175 // 176 //eliminates the island with a time spread 177 //higher than a given limit 178 // 179 //else if( fIslandCleaningMethod == kTiming ){ 180 else if( fIslandCleaningMethod == 0 ){ 181 while ((imgIsl=(MImgIsland*)Next())) { 182 183 //fIslCleanThreshold has different values, FIXME, put two variables 184 if(imgIsl->GetTimeSpread() > fIslCleanThres) 202 203 //////////////////////////////////////////////////// 204 // 205 // ISLAND SIZE OVER EVENT SIZE ISLAND CLEANING 206 // eliminates the island with an island size over event size 207 // lower than a given limit 208 /////////////////////////////////////////////////// 209 else if ( fIslandCleaningMethod == 2 ) { 210 Float_t size = 0; 211 while ((imgIsl=(MImgIsland*)Next())) { 212 size += imgIsl->GetSizeIsl(); 213 } 214 215 Next.Reset(); 216 while ((imgIsl=(MImgIsland*)Next())) { 217 218 if(imgIsl->GetSizeIsl()/size < fIslCleanThres) 185 219 { 186 220 pixNum = imgIsl->GetPixNum(); … … 193 227 } 194 228 } 195 } 196 197 } 198 199 // 200 //eliminates all the islands except the continent, 201 //i.e. the larger island in the event 202 // 229 } 230 } 231 232 233 //////////////////////////////////////////////////// 234 // 235 // CONTINENT ISLAND CLEANING 236 // eliminates all the islands except the continent 237 // i.e. the larger island in the event 238 // according the number of pixels 239 /////////////////////////////////////////////////// 203 240 else if( fIslandCleaningMethod == 3 ){ 204 241 Int_t i = 0; … … 219 256 } 220 257 258 259 //////////////////////////////////////////////////// 260 // 261 // LARGER and SECOND LARGER ISLAND CLEANING 262 // eliminates all the islands except the two biggest 263 // ones according size 264 // 265 /////////////////////////////////////////////////// 266 else if( fIslandCleaningMethod == 4 ){ 267 268 Int_t i = 0; 269 Int_t islnum = fIsl->GetIslNum(); 270 Float_t islSize[islnum]; 271 Int_t islIdx[islnum]; 272 273 for (Int_t j = 0; j<islnum ; j++){ 274 islSize[j] = -1; 275 islIdx[j] = -1; 276 } 277 278 while ((imgIsl=(MImgIsland*)Next())) { 279 280 islSize[i] = imgIsl->GetSizeIsl(); 281 i++; 282 } 283 284 285 TMath::Sort(islnum, islSize, islIdx, kTRUE); 286 287 i = 0; 288 Next.Reset(); 289 while ((imgIsl=(MImgIsland*)Next())) { 290 if (islnum > 1 && islIdx[0]!=i && islIdx[1]!=i){ 291 292 //cout << "removed " << i << " isl 0" << islIdx[0] << " isl 1" << islIdx[1] << endl; 293 294 pixNum = imgIsl->GetPixNum(); 295 296 for(Int_t k = 0; k<pixNum; k++) 297 { 298 idPix = imgIsl->GetPixList(k); 299 MCerPhotPix &pix = (*fEvt)[idPix]; 300 pix.SetPixelUnused(); 301 } 302 } 303 i++; 304 } 305 } 306 307 308 309 /////////////////////////////////////////////////////// 310 // 311 // LARGER and SECOND LARGER ISLAND CLEANING II 312 // eliminates all the islands except the two biggest 313 // ones according size, if the second one almost has 314 // the 80% of the size of the biggest one 315 // 316 // 317 ////////////////////////////////////////////////////// 318 else if( fIslandCleaningMethod == 5 ){ 319 320 Int_t i = 0; 321 Int_t islnum = fIsl->GetIslNum(); 322 Float_t islSize[islnum]; 323 Int_t islIdx[islnum]; 324 325 for (Int_t j = 0; j<islnum ; j++){ 326 islSize[j] = -1; 327 islIdx[j] = -1; 328 } 329 330 while ((imgIsl=(MImgIsland*)Next())) { 331 332 islSize[i] = imgIsl->GetSizeIsl(); 333 i++; 334 } 335 336 337 TMath::Sort(islnum, islSize, islIdx, kTRUE); 338 339 i = 0; 340 Next.Reset(); 341 while ((imgIsl=(MImgIsland*)Next())) { 342 343 if (islnum > 1 && islIdx[0]!=i && islIdx[1]!=i){ 344 345 pixNum = imgIsl->GetPixNum(); 346 347 for(Int_t k = 0; k<pixNum; k++) 348 { 349 idPix = imgIsl->GetPixList(k); 350 MCerPhotPix &pix = (*fEvt)[idPix]; 351 pix.SetPixelUnused(); 352 } 353 } 354 else if(islnum>1 && islSize[islIdx[1]]<0.6*islSize[islIdx[0]]){ 355 356 pixNum = imgIsl->GetPixNum(); 357 358 for(Int_t k = 0; k<pixNum; k++) 359 { 360 idPix = imgIsl->GetPixList(k); 361 MCerPhotPix &pix = (*fEvt)[idPix]; 362 pix.SetPixelUnused(); 363 } 364 } 365 i++; 366 } 367 } 368 221 369 fEvt->SetReadyToSave(); 222 370
Note:
See TracChangeset
for help on using the changeset viewer.