Ignore:
Timestamp:
08/02/04 09:27:43 (20 years ago)
Author:
aliu
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mifae
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtemp/mifae/Changelog

    r4431 r4439  
    1919                                                 -*-*- END OF LINE -*-*-
    2020
     21  2004/08/02 Ester Aliu Fusté
     22    * library/MIslands.[h, cc], MIslandCalc.[h,cc]
     23     - Added the variables: distance (Dist) between the larger island and
     24       the other ones, width and lenght of the larger island.
     25     
     26     - Written pointers instead of vectors
     27     
    2128  2004/07/28 Javi Lopez
    2229    * script/

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

    r4423 r4439  
    11/* ======================================================================== *\
     2!
     3! *
     4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
     5! * Software. It is distributed to you in the hope that it can be a useful
     6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
     7! * It is distributed WITHOUT ANY WARRANTY.
     8! *
     9! * Permission to use, copy, modify and distribute this software and its
     10! * documentation for any purpose is hereby granted without fee,
     11! * provided that the above copyright notice appear in all copies and
     12! * that both that copyright notice and this permission notice appear
     13! * in supporting documentation. It is provided "as is" without express
     14! * or implied warranty.
     15! *
    216!
    317!
    418!   Author(s): Ester Aliu, 2/2004 <aliu@ifae.es>
    5 !   
    6 !   Last Update: 6/2004
     19|
     20!   Last Update: 7/2004
     21!
     22!
     23!   Copyright: MAGIC Software Development, 2000-2004
     24!
    725!
    826\* ======================================================================== */
    … …  
    1533// of the events such as:
    1634//
    17 //   - fIslNum             //  number of islands
    18 //   - fIslId[577]         //  island Id
    19 //   - fPixNum[20]         //  number of pixels in the island
    20 //   - fSigToNoise[20]     //  signal to noise of the island
    21 //   - fTime[20][577]      //  mean of the arrival time 
    22 //   - fTimeSpread[20]     //  mean arrival time spread of the island
    23 //   - fMeanX[20]          //  mean X position of the island
    24 //   - fMeanY[20]          //  mean Y position of the island
    25 //   - fDist[20]           //  dist between an island and the continent
     35//   - fIslNum                 //  number of islands
     36//   - fIslId[577]             //  island Id
     37//   - fPixNum[numisl]         //  number of pixels in the island
     38//   - fSigToNoise[numisl]     //  signal to noise of the island
     39//   - ftime[numisl][577]      //  mean of the arrival time 
     40//   - fTimeSpread[numisl]     //  mean arrival time spread of the island
     41//   - fMeanX[numisl]          //  mean X position of the island
     42//   - fMeanY[numisl]          //  mean Y position of the island
     43//   - fDist[numisl]           //  dist between an island and the continent
     44//   - fLength                 //  major axis of the larger island ellipse
     45//   - fWidth                  //  minor axis of the larger island ellipse
     46//   - fDistL[numisl]          //  dist divided by lenght of the larger island
     47//   - fDistW[numisl]          //  dist divided by width of the larger island
    2648//
    2749// Input Containers:
    … …  
    124146
    125147
    126 Int_t MIslandCalc::Process()
    127 {
    128  
    129   if (fIslandAlgorithm == 1)
    130     Calc1();
    131  
    132   if (fIslandAlgorithm == 2)
    133     Calc2();
     148Int_t MIslandCalc::Process(){
     149 
     150  IslandPar();
     151 
    134152  return kTRUE; 
    135153   
    … …  
    137155
    138156
    139 Int_t MIslandCalc::Calc1(){
     157Int_t MIslandCalc::IslandPar(){
     158
     159  //calculates all the island parameters
     160
     161  const Int_t nPix=577;
     162  const Int_t nVect=20;
     163  Int_t numisl;
     164
     165  Int_t** vect;
     166  vect = new Int_t*[nVect];
     167  for(Int_t i=0;i<nVect;i++)
     168    vect[i]= new Int_t[nPix];
     169 
     170  Int_t num[nVect];
     171  // num = new Int_t*[nVect];
     172 
     173  if (fIslandAlgorithm == 1)
     174    Calc1(numisl,nVect,nPix,vect,num);
     175  if (fIslandAlgorithm == 2)
     176    Calc2(numisl,nVect,nPix,vect,num);
     177 
     178  //set the number of islands in one event
     179 
     180  fIsl->SetIslNum(numisl);
     181
     182 
     183  //examine each island...
     184
     185  Float_t  noise;
     186  Float_t  signal;
     187
     188  /*
     189  Float_t** ftime;
     190  ftime = new Float_t*[numisl];
     191
     192  Int_t pixNumIsl = 0;
     193
     194  for(Int_t i=0;i<numisl;i++)
     195    {
     196      pixNumIsl = num[i+1];
     197      ftime[i]= new Float_t[pixNumIsl];
     198    }
     199
     200  Int_t** fIslId;
     201  fIslId = new Int_t*[numisl];
     202  for(Int_t i=0;i<numisl;i++)
     203    {
     204      pixNumIsl = num[i+1];
     205      fIslId[i]= new Int_t[pixNumIsl];
     206    }
     207  */
     208 
     209  Int_t fPixNum[numisl];
     210  Float_t fSigToNoise[numisl];
     211  Float_t time[nPix];
     212  Float_t timeVariance[numisl];
     213  Float_t meanX[numisl];
     214  Float_t meanY[numisl];
     215  Float_t dist[numisl];
     216  Float_t distL[numisl];
     217  Float_t distW[numisl];
     218  Float_t size, sizeLargeIsl, length, width;
     219
     220
     221  //reset the "sets" functions
     222  if (numisl <1)
     223    fIsl->SetIslNum(0);
     224 
     225  for(Int_t i = 0; i<10 ;i++){
     226    for(Int_t idx = 0; idx<nPix; idx++)
     227      {
     228        if (i == 0)
     229          fIsl->SetIslId(idx, -1);
     230        fIsl->SetArrivalTime(i, idx, -1 );
     231      }
     232  }
     233 
     234  Float_t X = 0;
     235  Float_t Y = 0;
     236  sizeLargeIsl = 0;
     237
     238
     239  for(Int_t i = 1; i<=numisl ; i++)
     240    {
     241      Int_t n = 0;
     242      //Int_t ncore = 0;
     243     
     244      Float_t MIN = 10000.;
     245      Float_t MAX = 0.;
     246     
     247      signal = 0;
     248      noise = 0;
     249
     250      size = 0;
     251      meanX[i-1] = 0;
     252      meanY[i-1] = 0;
     253      dist[i-1] = 0;
     254
     255      fPixNum[i-1] = 0;
     256      timeVariance[i-1] = 0;
     257     
     258
     259      for(Int_t idx=0 ; idx<nPix ; idx++)
     260        {
     261         
     262          MCerPhotPix *pix = fEvt->GetPixById(idx);
     263          const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()];
     264          const MPedestalPix &ped  = (*fPed)[idx];
     265          const MArrivalTimePix &timepix = (*fTime)[idx];
     266          const Float_t nphot = pix->GetNumPhotons();
     267
     268          if (pix == NULL) break;
     269           
     270          if (vect[i][idx]==1){
     271           
     272            fPixNum[i-1]++;
     273            signal += nphot * (fCam->GetPixRatio(idx));
     274            noise += pow(ped.GetPedestalRms(),2);
     275           
     276            size += nphot;
     277            if (i == 1)
     278              sizeLargeIsl += nphot;
     279
     280            meanX[i-1] += nphot * gpix2.GetX();
     281            meanY[i-1] += nphot * gpix2.GetY();
     282           
     283            time[i-1] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
     284                       
     285            //      ftime[i-1][n] = time[i-1];
     286            //  fIslId[i-1][n] = idx;
     287
     288            //calculates the time spread only for core pixels 
     289            if (fEvt->IsPixelCore(idx)){
     290             
     291              if (time[i-1] > MAX)
     292                MAX = time[i-1];
     293              if (time[i-1] < MIN)
     294                MIN = time[i-1];
     295              //  ncore++;
     296            }
     297           
     298            fIsl->SetIslId(idx, i-1);
     299            fIsl->SetArrivalTime(i-1, idx, time[n]);
     300
     301            n++;
     302          }
     303         
     304        }
     305
     306      meanX[i-1] /= size;
     307      meanY[i-1] /= size;
     308 
     309
     310      if (i == 1){
     311        X = meanX[i-1];
     312        Y = meanY[i-1];
     313      }
     314 
     315      dist[i-1] = TMath::Power(meanX[i-1]-X,2) + TMath::Power(meanY[i-1]-Y,2);
     316      dist[i-1] = TMath::Sqrt(dist[i-1]);
     317     
     318      timeVariance[i-1] = MAX-MIN;
     319     
     320      fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
     321     
     322    }
     323
     324  //fIsl->SetIslId(fIslId);
     325  //fIsl->SetArrivalTime(ftime);
     326  fIsl->SetPixNum(fPixNum);
     327  fIsl->SetSigToNoise(fSigToNoise);
     328  fIsl->SetTimeSpread(timeVariance);
     329  fIsl->SetMeanX(meanX);
     330  fIsl->SetMeanY(meanY);
     331  fIsl->SetDist(dist);
     332 
     333
     334  //Length and Width of the larger island according the definition of the hillas parameters
     335 
     336  // calculate 2nd moments
     337  // ---------------------
     338  Double_t corrxx=0;                               // [m^2]
     339  Double_t corrxy=0;                               // [m^2]
     340  Double_t corryy=0;                               // [m^2]
     341 
     342  for(Int_t idx=0 ; idx<nPix ; idx++)
     343    {
     344      MCerPhotPix *pix = fEvt->GetPixById(idx);
     345      const MGeomPix &gpix3 = (*fCam)[pix->GetPixId()];
     346      const Float_t nphot = pix->GetNumPhotons();
     347     
     348      if (pix == NULL) break;
     349     
     350      if (vect[1][idx]==1){
     351       
     352        const Float_t dx = gpix3.GetX() - X;     // [mm]
     353        const Float_t dy = gpix3.GetY() - Y;     // [mm]
     354       
     355       
     356        corrxx += nphot * dx*dx;                     // [mm^2]
     357        corrxy += nphot * dx*dy;                     // [mm^2]
     358        corryy += nphot * dy*dy;                     // [mm^2]
     359       
     360      }   
     361    }
     362 
     363  // calculate the hillas parameters Width and Length
     364 
     365  const Double_t d0    = corryy - corrxx;
     366  const Double_t d1    = corrxy*2;
     367  const Double_t d2    = d0 + TMath::Sqrt(d0*d0 + d1*d1);
     368  const Double_t tand  = d2 / d1;
     369  const Double_t tand2 = tand*tand;
     370 
     371  const Double_t s2 = tand2+1;
     372 
     373  const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/sizeLargeIsl;
     374  const Double_t axis2 = (tand2*corrxx - d2 + corryy)/s2/sizeLargeIsl;
     375 
     376  //
     377  // fLength^2 is the second moment along the major axis of the ellipse
     378  // fWidth^2  is the second moment along the minor axis of the ellipse
     379  //
     380  // From the algorithm we get: fWidth <= fLength is always true
     381  //
     382  // very small numbers can get negative by rounding
     383  //
     384  length = axis1<0 ? 0 : TMath::Sqrt(axis1);  // [mm]
     385  width  = axis2<0 ? 0 : TMath::Sqrt(axis2);  // [mm]
     386 
     387  fIsl->SetLength(length);
     388  fIsl->SetWidth(width);
     389   
     390  // for(Int_t i = 1; i<=numisl ; i++){
     391
     392  // fIsl->SetDistL(fIsl->GetDist(i-1)/length, i-1);
     393  // fIsl->SetDistW(fIsl->GetDist(i-1)/width, i-1);
     394  // }
     395
     396  for(Int_t i = 1; i<=numisl ; i++){
     397   
     398    distL[i-1]=dist[i-1]/length;
     399    distW[i-1]=dist[i-1]/width;
     400   
     401  }
     402
     403  fIsl->SetDistL(distL);
     404  fIsl->SetDistW(distW);
     405
     406  fIsl->SetReadyToSave();
     407
     408  /*delete [] vect;
     409  delete [] num;
     410  delete [] ftime;
     411  delete [] fIslId;
     412  */
     413  return kTRUE; 
     414}
     415
     416//------------------------------------------------------------------------------------------
     417void MIslandCalc::Calc1(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){
    140418 
    141419 
    … …  
    148426  /////////////////////////////
    149427 
    150   Float_t  noise;
    151   Float_t  signal;
    152 
    153   Int_t    npix = 577;
    154 
    155428  Int_t    sflag;
    156429  Int_t    control;
    157430 
    158431  Int_t    nvect = 0;
    159   Int_t    fIslNum = 0;
    160  
    161   Int_t    vect[50][577];
    162 
    163   Int_t    zeros[50];
    164  
    165   for(Int_t m = 0; m < 50 ; m++)
     432 
     433  numisl = 0;
     434
     435  Int_t    zeros[nv];
     436 
     437  for(Int_t m = 0; m < nv ; m++)
    166438    for(Int_t n = 0; n < npix ; n++)
    167439        vect[m][n] = 0;
    168440   
    169   for(Int_t n = 0; n < 50 ; n++)
     441  for(Int_t n = 0; n < nv ; n++)
    170442    zeros[n] = 0;
    171443 
    … …  
    213485    }
    214486 
    215   fIslNum = nvect;
     487  numisl = nvect;
    216488 
    217489 
    … …  
    241513                  zeros[j] = 1;
    242514                }       
    243               fIslNum = fIslNum-1;         
     515              numisl = numisl-1;           
    244516            }
    245517         
    … …  
    265537              if (vect[l][k] == 1)
    266538                numpixels++;
     539                       
    267540            }
    268541          if (numpixels>pixMAX)
    … …  
    273546          l++;
    274547        }
     548      num[i] = numpixels;
     549
    275550    }
    276551 
    277552  //the larger island will correspond to the 1st component of the vector
     553
     554  num[nvect +1] = num[1];
     555  num[1] = num[idMAX];
     556  num[idMAX]=num[1];
    278557
    279558  for(Int_t k = 0; k<npix; k++)
    … …  
    283562      vect[idMAX][k] = vect[nvect+1][k];
    284563    }
    285 
    286 
    287   //set the number of islands in one event
    288   fIsl->SetIslNum(fIslNum);
    289 
    290   //examine each island...
    291   Int_t fPixNum[fIslNum];
    292   Float_t fSigToNoise[fIslNum];
    293   Float_t time[577];
    294   Float_t timeVariance[fIslNum];
    295   Float_t size, meanX, meanY, dist;
    296 
    297   //reset the "sets" functions
    298   if (fIslNum <1)
    299     fIsl->SetIslNum(0);
    300 
    301   for(Int_t i = 0; i<20 ;i++)
    302     {
    303       fIsl->SetPixNum(i,-1);
    304       fIsl->SetSigToNoise(i,-1);
    305       fIsl->SetTimeSpread(i,-1);
    306       fIsl->SetMeanX(i-1, -10000);
    307       fIsl->SetMeanY(i-1, -10000);
    308       fIsl->SetDist(i-1, -1);
    309 
    310       for(Int_t idx = 0; idx<npix; idx++)
    311         {
    312           fIsl->SetIslId(idx, -1);
    313           fIsl->SetArrivalTime(i, idx, -1 );
    314         }
    315     }
    316  
    317 
    318   Float_t X = 0;
    319   Float_t Y = 0;
    320    
    321    for(Int_t i = 1; i<=fIslNum ; i++)
    322     {
    323       Int_t n = 0;
    324       Int_t ncore = 0;
    325      
    326       Float_t MIN = 10000;
    327       Float_t MAX = 0;
    328 
    329       signal = 0;
    330       noise = 0;
    331 
    332       size = 0;
    333       meanX = 0;
    334       meanY = 0;
    335       dist = 0;
    336 
    337       fPixNum[i-1] = 0;
    338       timeVariance[i-1] = 0;
    339      
    340       for(Int_t idx=0 ; idx<npix ; idx++)
    341         {
    342          
    343           MCerPhotPix *pix = fEvt->GetPixById(idx);
    344           const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()];
    345           const MPedestalPix &ped  = (*fPed)[idx];
    346           const MArrivalTimePix &timepix = (*fTime)[idx];
    347           const Float_t nphot = pix->GetNumPhotons();
    348 
    349           if (pix == NULL) break;
    350            
    351           if (vect[i][idx]==1){
    352            
    353             fPixNum[i-1]++;
    354             signal += nphot * (fCam->GetPixRatio(idx));
    355             noise += pow(ped.GetPedestalRms(),2);
    356            
    357             size += nphot;
    358             meanX += nphot * gpix2.GetX();
    359             meanY += nphot * gpix2.GetY();
    360            
    361             time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
    362            
    363            
    364             if (fEvt->IsPixelCore(idx)){
    365              
    366               if (time[n] > MAX)
    367                 MAX = time[n];
    368               if (time[n] < MIN)
    369                 MIN = time[n];
    370              
    371               ncore++;
    372             }
    373            
    374             fIsl->SetIslId(idx, i-1);
    375             fIsl->SetArrivalTime(i-1, idx, time[n]);
    376                
    377             n++;
    378           }
    379          
    380         }
    381 
    382       meanX /= size;
    383       meanY /= size;
    384  
    385       if (i ==1){
    386         X = meanX;
    387         Y = meanY;
    388       }
    389    
    390       dist = TMath::Power(meanX-X,2) + TMath::Power(meanY-Y,2);
    391       dist = TMath::Sqrt(dist);
    392 
    393       timeVariance[i-1] = (MAX - MIN);
    394 
    395       fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
    396 
    397       fIsl->SetPixNum(i-1,fPixNum[i-1]);
    398       fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
    399       fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
    400       fIsl->SetMeanX(i-1, meanX);
    401       fIsl->SetMeanY(i-1, meanY);
    402       fIsl->SetDist(i-1, dist);
    403      
    404     }
    405  
    406   fIsl->SetReadyToSave();
    407  
    408   return kTRUE; 
    409564}
    410565
    411 
    412 Int_t MIslandCalc::Calc2(){
    413  
     566//------------------------------------------------------------------------------------------
     567
     568void MIslandCalc::Calc2(Int_t& numisl, const Int_t nv, const Int_t npix, Int_t** vect, Int_t* num){ 
     569
    414570 
    415571  /////////////////////////////
    … …  
    421577  /////////////////////////////
    422578 
    423   Float_t  noise;
    424   Float_t  signal;
    425 
    426   Int_t    npix = 577;
    427 
    428579  Int_t    sflag;
    429580  Int_t    control;
    430581 
    431582  Int_t    nvect = 0;
    432   Int_t    fIslNum = 0;
     583  numisl = 0;
    433584 
    434   Int_t    vect[50][577];
    435 
    436   Int_t    zeros[50];
    437  
    438   Int_t kk[577];
    439 
    440   for(Int_t m = 0; m < 50 ; m++)
     585  Int_t    zeros[nv];
     586 
     587  Int_t kk[npix];
     588
     589  for(Int_t m = 0; m < nv ; m++)
    441590    for(Int_t n = 0; n < npix ; n++)
    442591        vect[m][n] = 0;
    443592   
    444   for(Int_t n = 0; n < 50 ; n++)
     593  for(Int_t n = 0; n < nv ; n++)
    445594    zeros[n] = 0;
    446595 
    447   for(Int_t n = 0; n < 577 ; n++)
     596  for(Int_t n = 0; n < npix ; n++)
    448597    kk[n] = 0;
    449598 
    … …  
    511660    }
    512661 
    513   fIslNum = nvect;
     662  numisl = nvect;
    514663 
    515664 
    … …  
    539688                  zeros[j] = 1;
    540689                }       
    541               fIslNum = fIslNum-1;         
     690              numisl = numisl-1;           
    542691            }
    543692         
    … …  
    570719          l++;
    571720        }
     721      num[i] = numpixels;
    572722    }
    573723 
    574724 
    575725  //the larger island will correspond to the 1st component of the vector
     726
     727  num[nvect +1] = num[1];
     728  num[1] = num[idMAX];
     729  num[idMAX]=num[1];
    576730
    577731  for(Int_t k = 0; k<npix; k++)
    … …  
    581735      vect[idMAX][k] = vect[nvect+1][k];
    582736    }
    583  
    584 
    585   //set the number of islands in one event
    586   fIsl->SetIslNum(fIslNum);
    587 
    588   //examine each island...
    589   Int_t fPixNum[fIslNum];
    590   Float_t fSigToNoise[fIslNum];
    591   Float_t time[577];
    592   Float_t timeVariance[fIslNum];
    593   Float_t size, meanX, meanY, dist;
    594 
    595   //reset the "sets" functions
    596   if (fIslNum <1)
    597     fIsl->SetIslNum(0);
    598 
    599   for(Int_t i = 0; i<20 ;i++)
    600     {
    601       fIsl->SetPixNum(i,-1);
    602       fIsl->SetSigToNoise(i,-1);
    603       fIsl->SetTimeSpread(i,-1);
    604       fIsl->SetMeanX(i-1, -10000);
    605       fIsl->SetMeanY(i-1, -10000);
    606       fIsl->SetDist(i-1, -1);
    607  
    608       for(Int_t idx = 0; idx<npix; idx++)
    609         {
    610           fIsl->SetIslId(idx, -1);
    611           fIsl->SetArrivalTime(i, idx, -1 );
    612         }
    613     }
    614  
    615 
    616   Float_t X = 0;
    617   Float_t Y = 0;
    618    
    619    for(Int_t i = 1; i<=fIslNum ; i++)
    620     {
    621       Int_t n = 0;
    622       Int_t ncore = 0;
    623      
    624       Float_t MIN = 10000;
    625       Float_t MAX = 0;
    626 
    627       signal = 0;
    628       noise = 0;
    629 
    630       size = 0;
    631       meanX = 0;
    632       meanY = 0;
    633       dist = 0;
    634 
    635       fPixNum[i-1] = 0;
    636       timeVariance[i-1] = 0;
    637      
    638       for(Int_t idx=0 ; idx<npix ; idx++)
    639         {
    640          
    641           MCerPhotPix *pix = fEvt->GetPixById(idx);
    642           const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()];
    643           const MPedestalPix &ped  = (*fPed)[idx];
    644           const MArrivalTimePix &timepix = (*fTime)[idx];
    645           const Float_t nphot = pix->GetNumPhotons();
    646 
    647           if (pix == NULL) break;
    648            
    649           if (vect[i][idx]==1){
    650            
    651             fPixNum[i-1]++;
    652             signal += nphot * (fCam->GetPixRatio(idx));
    653             noise += pow(ped.GetPedestalRms(),2);
    654            
    655             size += nphot;
    656             meanX += nphot * gpix2.GetX();
    657             meanY += nphot * gpix2.GetY();
    658            
    659             time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
    660            
    661            
    662             if (fEvt->IsPixelCore(idx)){
    663              
    664               if (time[n] > MAX)
    665                 MAX = time[n];
    666               if (time[n] < MIN)
    667                 MIN = time[n];
    668              
    669               ncore++;
    670             }
    671            
    672             fIsl->SetIslId(idx, i-1);
    673             fIsl->SetArrivalTime(i-1, idx, time[n]);
    674                
    675             n++;
    676           }
    677          
    678         }
    679 
    680       meanX /= size;
    681       meanY /= size;
    682  
    683       if (i ==1){
    684         X = meanX;
    685         Y = meanY;
    686       }
    687    
    688       dist = TMath::Power(meanX-X,2) + TMath::Power(meanY-Y,2);
    689       dist = TMath::Sqrt(dist);
    690 
    691       timeVariance[i-1] = (MAX - MIN);
    692 
    693       fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
    694 
    695       fIsl->SetPixNum(i-1,fPixNum[i-1]);
    696       fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
    697       fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
    698       fIsl->SetMeanX(i-1, meanX);
    699       fIsl->SetMeanY(i-1, meanY);
    700       fIsl->SetDist(i-1, dist);
    701      
    702     }
    703 
    704      
    705   fIsl->SetReadyToSave();
    706  
    707   return 1; 
     737
    708738}
     739

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

    r4286 r4439  
    3636    Int_t PreProcess(MParList *plist);
    3737    Int_t Process();
    38     Int_t Calc1();                  // algorithm of counting islands #1
    39     Int_t Calc2();                  // algorithm of counting islands #2   
     38    Int_t IslandPar();               //
     39    void  Calc1(Int_t&,const Int_t,const Int_t,Int_t**,Int_t*);    // algorithm of counting islands #1
     40    void  Calc2(Int_t&,const Int_t,const Int_t,Int_t**,Int_t*);    // algorithm of counting islands #2
     41   
    4042           
    4143 public:

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

    r4418 r4439  
    3030MIslands::~MIslands()
    3131{
     32 
    3233}
     34
     35
     36// --------------------------------------------------------------------------
     37//
     38// Getter functions
     39//
     40
     41
     42// --------------------------------------------------------------------------
     43//
     44// Setter functions
     45//
     46
     47
     48
    3349
    3450// --------------------------------------------------------------------------
    … …  
    3652// Print the island parameters to *fLog
    3753//
    38 void MIslands::Print(Option_t *) const
     54void MIslands::Print(Option_t *opt=NULL) const
    3955{
    4056    *fLog << all;
    … …  
    4662        *fLog << "    + Pixel Number = " << fPixNum[i] << endl;
    4763        *fLog << "    + SigToNoise = " << fSigToNoise[i] << endl;
    48         *fLog << "    + TimeSpread = " << fTimeSpread[i] << endl;
    49         *fLog << "    + MeanX = " << fMeanX[i] << endl;
    50         *fLog << "    + MeanY = " << fMeanY[i] << endl;
    51         *fLog << "    + Dist = " << fDist[i] << endl;
     64        *fLog << "    + TimeSpread   [time slices]= " << fTimeSpread[i] << endl;
     65        *fLog << "    + MeanX   [mm]= " << fMeanX[i] << endl;
     66        *fLog << "    + MeanY   [mm]= " << fMeanY[i] << endl;
     67        *fLog << "    + Dist    [mm]= " << fDist[i] << endl;
     68        *fLog << "    + Length of the larger island  [mm]  = "  << fLength <<endl;
     69        *fLog << "    + Width of the larger island  [mm]  = "  << fWidth <<endl;
     70        *fLog << "    + DistL  = "  << fDistL[i] <<endl;
     71        *fLog << "    + DistW  = " << fDistW[i]  << endl;
    5272      }
    5373}

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

    r4418 r4439  
    1414#endif
    1515
     16#ifndef ROOT_TObjArray
     17#include <TObjArray.h>
     18#endif
    1619
    1720class MIslands : public MParContainer
    … …  
    1922private:
    2023    // for description see MIslands.cc
    21     Int_t fIslNum;                 //  number of islands
    22     Int_t fIslId[577];             //  island Id
    23     Int_t fPixNum[20];             //  number of pixels in the island
    24     Float_t fSigToNoise[20];       //  signal to noise of the island
    25     Float_t fTime[20][577];        //  mean of the arrival time 
    26     Float_t fTimeSpread[20];       //  mean arrival time spread of the island 
    27     Float_t fMeanX[20];            //  mean X position of the island
    28     Float_t fMeanY[20];            //  mean Y position of the island
    29     Float_t fDist[20];    //  dist between islands and continent(larger island)
    30 
     24    Int_t fIslNum;             //  number of islands
     25    //    Int_t** fIslId;      //[fIslNum]  island Id
     26    // TObjArray fIslId;       //  island Id
     27    Int_t fIslId[577];         //island Id
     28    Int_t* fPixNum;            //[fIslNum]  number of pixels in the island
     29    Float_t* fSigToNoise;      //[fIslNum]  signal to noise of the island
     30    //  Float_t** fTime;       //[fIslNum]  mean of the arrival time
     31    Float_t fTime[10][577];        //  mean of the arrival time
     32    Float_t* fTimeSpread;      //[fIslNum]  mean arrival time spread of the core pixels of the island 
     33    Float_t* fMeanX;           //[fIslNum]  mean X position of the island
     34    Float_t* fMeanY;           //[fIslNum]  mean Y position of the island
     35    Float_t* fDist;            //[fIslNum]  dist between islands and continent(larger island)
     36    Float_t fLength;           //  major axis of the larger island ellipse
     37    Float_t fWidth;            //  minor axis of the larger island ellipse
     38    Float_t* fDistL;           //[fIslNum] Dist of the island divided by Length of the larger island
     39    Float_t* fDistW;           //[fIslNum] Dist of the island divided by Width of the larger island
     40   
    3141public:
    3242    MIslands(const char *name=NULL, const char *title=NULL);
    3343    ~MIslands();
    3444
     45    // void Clear();
    3546    void Print(Option_t *opt=NULL) const;
     47   
     48    //getter methods
     49    Int_t    GetIslNum() const               { return fIslNum; }
     50    Int_t    GetIslId(Int_t idx)  { return fIslId[idx]; }
     51    //Int_t    GetIslId(Int_t isl, Int_t idx)  { return fIslId[isl][idx]; }
     52    // TObjArray GetIslId()                    {return fIslId;}
     53    Float_t  GetArrivalTime(Int_t isl, Int_t idx) { return fTime[isl][idx]; }     
     54    //TObjArray GetArrivalTime()               { return fTime; }     
     55    Int_t    GetPixNum(Int_t isl)            { return fPixNum[isl]; }
     56    Float_t  GetSigToNoise(Int_t isl)        { return fSigToNoise[isl]; }     
     57    Float_t  GetTimeSpread(Int_t isl)        { return fTimeSpread[isl];}         
     58    Float_t  GetMeanX(Int_t isl)             { return fMeanX[isl];}
     59    Float_t  GetMeanY(Int_t isl)             { return fMeanY[isl];}
     60    Float_t  GetDist(Int_t isl)              { return fDist[isl]; }
     61    Float_t  GetDistL(Int_t isl)             { return fDistL[isl]; }
     62    Float_t  GetDistW(Int_t isl)             { return fDistW[isl]; }
     63   
     64    Float_t  GetLength() const               { return fLength; }
     65    Float_t  GetWidth() const                { return fWidth; }
    3666
    37     Int_t    GetIslNum() const                        { return fIslNum; }
    38     Int_t    GetIslId(Int_t idx) const                { return fIslId[idx]; }
    39     Int_t    GetPixNum(Int_t isl)                     { return fPixNum[isl]; }
    40     Float_t  GetSigToNoise(Int_t isl)                 { return fSigToNoise[isl]; }
    41     Float_t  GetArrivalTime(Int_t isl, Int_t idx) { return fTime[isl][idx]; }           
    42     Float_t  GetTimeSpread(Int_t isl)             { return fTimeSpread[isl];}         
    43     Float_t  GetMeanX(Int_t isl)                  { return fMeanX[isl];}
    44     Float_t  GetMeanY(Int_t isl)                  { return fMeanY[isl];}
    45     Float_t  GetDist(Int_t isl)                   { return fDist[isl]; }
     67    //setter functions   
     68    void     SetIslNum(Int_t nisl)           { fIslNum = nisl; }
     69    void     SetIslId(Int_t idx, Int_t isl)           { fIslId[idx] = isl; }
     70   
     71    // void     SetIslId(Int_t** vect)           { fIslId = vect; }
     72    void     SetArrivalTime(Int_t isl, Int_t idx, Float_t val)   { fTime[isl][idx] = val;}
     73    //  void     SetArrivalTime(Float_t** vect)   { fTime = vect;}
     74    void     SetPixNum(Int_t* npix)          { fPixNum = npix; }
     75    void     SetSigToNoise(Float_t* val)     { fSigToNoise = val; }
     76    void     SetTimeSpread(Float_t* val)     { fTimeSpread = val; }
     77    void     SetMeanX( Float_t* val)         { fMeanX = val; }
     78    void     SetMeanY(Float_t* val)          { fMeanY = val; }
     79    void     SetDist(Float_t* val)           { fDist = val; }
     80    void     SetDistL(Float_t* val)          { fDistL=val; }
     81    void     SetDistW(Float_t* val)          { fDistW=val; }
    4682
    47     void     SetIslNum(Int_t nisl)                    { fIslNum = nisl; }
    48     void     SetIslId(Int_t idx, Int_t isl)           { fIslId[idx] = isl; }
    49     void     SetPixNum(Int_t isl, Int_t npix)         { fPixNum[isl] = npix; }
    50     void     SetSigToNoise(Int_t isl, Float_t val)    { fSigToNoise[isl] = val; }
    51     void     SetArrivalTime(Int_t isl, Int_t idx, Float_t val)   { fTime[isl][idx] = val;}
     83    void     SetLength(Float_t val)          { fLength=val; }
     84    void     SetWidth(Float_t val)           { fWidth=val; }
    5285
    53     void     SetTimeSpread(Int_t isl, Float_t val) { fTimeSpread[isl] = val; }
    54     void     SetMeanX(Int_t isl, Float_t val) { fMeanX[isl]=val; }
    55     void     SetMeanY(Int_t isl, Float_t val) { fMeanY[isl]=val; }
    56     void     SetDist(Int_t isl, Float_t val) { fDist[isl]=val; }
     86   
     87   
    5788
    58 
    59     ClassDef(MIslands, 1) // Storage Container for Island Parameters
     89    ClassDef(MIslands, 2) // Storage Container for Island Parameters
    6090};
    6191
Note: See TracChangeset for help on using the changeset viewer.