Ignore:
Timestamp:
11/24/04 16:33:39 (20 years ago)
Author:
aliu
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mifae
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mtemp/mifae/Changelog

    r5379 r5466  
    1818
    1919                                                 -*-*- 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
    2028
    2129 2004/11/11  Ester Aliu
  • trunk/MagicSoft/Mars/mtemp/mifae/library/MImgIsland.cc

    r5379 r5466  
    3737    fName  = name  ? name  : "MImgIsland";
    3838    fTitle = title ? title : "Storage container for one island information";
    39      
     39    //    cout << "Creo 1 isla" << endl;     
    4040    Reset();
     41}
     42MImgIsland::~MImgIsland()
     43{
     44  //  cout << "Destruyo 1 isla" << endl;
    4145}
    4246
     
    4448{
    4549  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;
    5363 
    5464  fPixList.Reset(-1);
    5565  fPeakPulse.Reset(-1);
     66   
    5667}
    5768
  • trunk/MagicSoft/Mars/mtemp/mifae/library/MIslandsCalc.cc

    r5379 r5466  
    3535//   - fPixNum                 //  number of pixels in the island
    3636//   - fSigToNoise             //  signal to noise of the island
    37 //   - fTimeSpread             //  mean arrival time spread of the island
     37//   - fTimeSpread             //  "time" of the island
    3838//   - fMeanX                  //  mean X position of the island
    3939//   - fMeanY                  //  mean Y position of the island
    4040//   - fDist                   //  dist between an island and the continent
    41 //   - fLength                 //  major axis of the larger island ellipse
    42 //   - fWidth                  //  minor axis of the larger island ellipse
     41//   - fLength                 //  major axis of the island ellipse
     42//   - fWidth                  //  minor axis of the island ellipse
    4343//   - fDistL                  //  dist divided by lenght of the larger island
    4444//   - fDistW                  //  dist divided by width of the larger island
     
    192192  Float_t meanX[numisl];
    193193  Float_t meanY[numisl];
     194  Float_t length[numisl];
     195  Float_t width[numisl];
    194196  Float_t dist[numisl];
    195197  Float_t distL[numisl];
     
    197199  Float_t distS[numisl];
    198200
    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
    203204  sizeLargeIsl = 0;
    204205  alphaW = 0;
     
    211212     
    212213      imgIsl->InitSize(num[i]);
    213 
     214     
    214215      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
    219226      signal = 0;
    220227      noise = 0;
    221 
     228     
    222229      size[i-1] = 0;
    223230      meanX[i-1] = 0;
    224231      meanY[i-1] = 0;
    225232      dist[i-1] = 0;
    226 
     233     
    227234      PixelNumIsl[i-1] = 0;
    228235      timeVariance[i-1] = 0;
    229      
    230      
     236           
    231237      for(Int_t idx=0 ; idx<nPix ; idx++)
    232238        {
     
    247253            if (i == 1)
    248254              sizeLargeIsl += nphot;
    249 
     255           
    250256            meanX[i-1] += nphot * gpix2.GetX();
    251257            meanY[i-1] += nphot * gpix2.GetY();
    252258           
    253259            time[i-1] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
    254            
     260            
    255261            imgIsl->SetPixList(PixelNumIsl[i-1]-1,pix->GetPixId());
    256262            imgIsl->SetPeakPulse(PixelNumIsl[i-1]-1,time[i-1]);                 
     
    259265            if (fEvt->IsPixelCore(idx)){
    260266             
    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            }
    268279            n++;
    269280          }       
    270281        } 
    271                    
     282     
    272283      meanX[i-1] /= size[i-1];
    273284      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);
    282287      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;
    285295     
    286296      //noise = 0, in the case of MC w/o noise
    287297      if (noise == 0) noise = 1;
    288 
     298     
    289299      SigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
    290 
     300     
    291301      imgIsl->SetPixNum(PixelNumIsl[i-1]);
    292302      imgIsl->SetSigToNoise(SigToNoise[i-1]);
     
    296306      imgIsl->SetDist(dist[i-1]);
    297307      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      }
    298450
    299451      fIsl->GetList()->Add(imgIsl);
    300452     
    301453    }
    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 
    441455  fIsl->SetAlphaW(alphaW/sizetot);   
    442456  //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;
    448463
    449464  return kTRUE; 
     
    481496  MCerPhotPix *pix;
    482497
    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
    486504  while ((pix=static_cast<MCerPhotPix*>(Next())))
    487505    {
     
    493511      if( fEvt->IsPixelUsed(idx))
    494512        {
    495           //cout << idx <<endl;
     513          //cout <<idx <<endl;
    496514          sflag = 0;
    497515         
     
    650668  MCerPhotPix *pix;
    651669
    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);
    654675 
    655676  while ((pix=static_cast<MCerPhotPix*>(Next0())))
     
    671692    }
    672693     
     694 
    673695  //2nd loop over all pixels
    674   MCerPhotEvtIter Next(fEvt, kFALSE);
    675  
     696  TIter Next(*fEvt);
     697   
    676698  while ((pix=static_cast<MCerPhotPix*>(Next())))
    677699    {
     
    713735  numisl = nvect;
    714736 
    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  }
    749781 
    750782  Int_t l = 1;
Note: See TracChangeset for help on using the changeset viewer.