Ignore:
Timestamp:
11/11/04 09:55:25 (20 years ago)
Author:
aliu
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mifae/library
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtemp/mifae/library/MControlPlots.cc

    r5170 r5379  
    1616!
    1717!   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)
    1920!   Copyright: MAGIC Software Development, 2000-2004
    2021!

  • trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc

    r5170 r5379  
    5454  fPixList.Reset(-1);
    5555  fPeakPulse.Reset(-1);
    56 
    5756}
    5857
    … …  
    105104  *fLog << inf << "  DistW = " << fDistW << endl;
    106105  *fLog << inf << "  DistS = " << fDistS << endl;
     106  *fLog << inf << "  Alpha = " << fAlpha << endl;
    107107
    108108}

  • trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.h

    r5170 r5379  
    2222  Float_t fSigToNoise;
    2323  Float_t fTimeSpread;
    24   Float_t fMeanX;
    25   Float_t fMeanY;
    26   Float_t fWidth;
    27   Float_t fLength;
    2824  Float_t fDist;
    2925  Float_t fDistL;
    3026  Float_t fDistW;
    3127  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
    3237
    3338  TArrayI fPixList;
    … …  
    4247  Float_t GetSigToNoise()  { return fSigToNoise; }
    4348  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; }
    4456  Float_t GetMeanX()       { return fMeanX; }
    4557  Float_t GetMeanY()       { return fMeanY; }
    4658  Float_t GetWidth()       { return fWidth; }
    4759  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; }
    5260
     61  // hillas src parameters
     62  Float_t GetAlpha()       { return fAlpha; }
     63 
    5364  void    InitSize(Int_t i);
    5465  UInt_t  GetSize() const { return fPixList.GetSize(); }
    … …  
    5970  void Reset();
    6071
    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; }
    7289 
    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);
    7592
    7693  //  void Paint(Option_t *opt=NULL);
    7794  void Print(Option_t *opt=NULL) const; 
    7895
    79   ClassDef(MImgIsland, 1) // Container that holds the island information
     96  ClassDef(MImgIsland, 2) // Container that holds the island information
    8097
    8198};

  • trunk/MagicSoft/Mars/mtemp/mifae/library/MIslands.h

    r5186 r5379  
    2323 
    2424  Int_t fIslNum;
     25  Float_t fAlphaW;
    2526
    2627 public:
    … …  
    3435  TList* GetList() const { return fIslands; }
    3536  UInt_t GetIslNum() const { return fIslands->GetSize(); } //number of islands
     37  Float_t GetAlphaW() const { return fAlphaW; }     //alpha weighted
    3638
    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; }
    3941  //  MHCamera& GetDisplay() { return fDisplay; }
    4042

  • trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc

    r5170 r5379  
    180180  fIsl->SetIslNum(numisl);
    181181
     182  //cout << "code numisl " << numisl << endl;
    182183  //examine each island...
    183184
    … …  
    196197  Float_t distS[numisl];
    197198
    198   Float_t size[numisl], sizeLargeIsl, length, width;
     199  Float_t size[numisl], sizeLargeIsl, length, width, distance, alpha, alphaW, sizetot;
    199200
    200201  Float_t X = 0;
    201202  Float_t Y = 0;
    202203  sizeLargeIsl = 0;
     204  alphaW = 0;
     205  sizetot = 0;
    203206
    204207  for(Int_t i = 1; i<=numisl ; i++)
    … …  
    292295      imgIsl->SetMeanY(meanY[i-1]);
    293296      imgIsl->SetDist(dist[i-1]);
     297      imgIsl->SetSizeIsl(size[i-1]);
    294298
    295299      fIsl->GetList()->Add(imgIsl);
    … …  
    333337  MImgIsland *imgIsl = new MImgIsland;
    334338  TIter Next(fIsl->GetList());
    335  
     339  //Next.Reset();
     340
    336341  Int_t i = 1;
    337342  while ((imgIsl=(MImgIsland*)Next())) {
    … …  
    344349   
    345350    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
    347356    const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/size[i-1];
    348357    const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/size[i-1];
    … …  
    359368    width  = axis2<0 ? 0 : TMath::Sqrt(axis2);  // [mm]
    360369   
    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   
    362425    distL[i-1]=dist[i-1]/length;
    363426    distW[i-1]=dist[i-1]/width;
    364427    distS[i-1]= dist[i-1]/size[i-1];
     428
    365429    imgIsl->SetLength(length);
    366430    imgIsl->SetWidth(width);
     431
    367432    imgIsl->SetDistL(distL[i-1]);
    368433    imgIsl->SetDistW(distW[i-1]);
    369434    imgIsl->SetDistS(distS[i-1]);
     435
     436    imgIsl->SetAlpha(alpha);
    370437    i++;
    371438  }
    372    
     439
     440
     441  fIsl->SetAlphaW(alphaW/sizetot);   
    373442  //fIsl->SetReadyToSave();
    374443
    … …  
    394463 
    395464  Int_t    sflag;
    396   Int_t    control;
     465  Int_t    control = 0;
    397466 
    398467  Int_t    nvect = 0;
    … …  
    455524  numisl = nvect;
    456525 
    457  
     526   
    458527  // Repeated Chain Corrections
    459528
    460529  Int_t jmin = 0;
    461530 
    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  }
    513571 
    514572  Int_t pixMAX = 0;
    … …  
    540598        }
    541599    }
    542  
     600
     601
    543602  //the larger island will correspond to the 1st component of the vector
    544603 
    … …  
    657716  // Repeated Chain Corrections
    658717 
    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  }
    715749 
    716750  Int_t l = 1;

  • trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsClean.cc

    r5170 r5379  
    144144Int_t MIslandsClean::Process()
    145145{
    146   //
    147   //eliminates the island with a signal-to-noise
    148   //lower than a given limit
    149   //
    150   //if ( fIslandCleaningMethod == kNoTiming ){
    151146 
    152147  MImgIsland *imgIsl = new MImgIsland;
    … …  
    156151  Int_t idPix = -1;
    157152
    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 ) {
    159186    while ((imgIsl=(MImgIsland*)Next())) {
    160187     
    … …  
    173200  } 
    174201 
    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)
    185219        {
    186220          pixNum = imgIsl->GetPixNum();
    … …  
    193227            }
    194228        }
    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  ///////////////////////////////////////////////////
    203240  else if( fIslandCleaningMethod == 3 ){
    204241    Int_t i = 0;
    … …  
    219256  }
    220257 
     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 
    221369  fEvt->SetReadyToSave();
    222370 
Note: See TracChangeset for help on using the changeset viewer.