Ignore:
Timestamp:
07/22/04 19:30:31 (20 years ago)
Author:
aliu
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mifae/library
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r4286 r4418  
    1717//   - fIslNum             //  number of islands
    1818//   - fIslId[577]         //  island Id
    19 //   - fPixNum[20];        //  number of pixels in the island
    20 //   - fSigToNoise[20];    //  signal to noise of the island
     19//   - fPixNum[20]         //  number of pixels in the island
     20//   - fSigToNoise[20]     //  signal to noise of the island
    2121//   - fTime[20][577]      //  mean of the arrival time 
    22 //   - fTimeSpread[20];    //  mean arrival time spread of the island
    23 //
     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
    2426//
    2527// Input Containers:
    … …  
    141143  //
    142144  //        ALGORITHM # 1
    143   // counts the number of algorithms as you can see in
     145  // counts the number of islands as you can see in
    144146  // the event display after doing the std image cleaning
    145147  //
    … …  
    215217 
    216218  // Repeated Chain Corrections
     219
    217220 
    218221  for(Int_t i = 1; i <= nvect; i++)
    … …  
    223226          for(Int_t k = 0; k < npix; k++)
    224227            {
    225              
    226228              if (vect[i][k] == 1 && vect[j][k] == 1)
    227229                {
    … …  
    246248 
    247249 
     250
    248251  Int_t l = 1;
    249  
     252  Int_t numpixels;
     253  Int_t pixMAX = 0;
     254  Int_t idMAX = 1;
     255
    250256  for(Int_t i = 1;  i<= nvect ; i++)
    251257    {
     258      numpixels = 0;
     259
    252260      if (zeros[i] == 0)
    253261        {
    … …  
    255263            {
    256264              vect[l][k] = vect[i][k];
     265              if (vect[l][k] == 1)
     266                numpixels++;
     267            }
     268          if (numpixels>pixMAX)
     269            {
     270              pixMAX = numpixels;
     271              idMAX = l;
    257272            }
    258273          l++;
    … …  
    260275    }
    261276 
    262  
     277  //the larger island will correspond to the 1st component of the vector
     278
     279  for(Int_t k = 0; k<npix; k++)
     280    {
     281      vect[nvect+1][k] = vect[1][k];
     282      vect[1][k] = vect[idMAX][k];
     283      vect[idMAX][k] = vect[nvect+1][k];
     284    }
     285
     286
    263287  //set the number of islands in one event
    264288  fIsl->SetIslNum(fIslNum);
    … …  
    269293  Float_t time[577];
    270294  Float_t timeVariance[fIslNum];
    271  
     295  Float_t size, meanX, meanY, dist;
     296
    272297  //reset the "sets" functions
    273298  if (fIslNum <1)
    … …  
    287312    }
    288313 
     314
     315  Float_t X = 0;
     316  Float_t Y = 0;
    289317   
    290318   for(Int_t i = 1; i<=fIslNum ; i++)
    … …  
    298326      signal = 0;
    299327      noise = 0;
     328
     329      size = 0;
     330      meanX = 0;
     331      meanY = 0;
     332      dist = 0;
     333
    300334      fPixNum[i-1] = 0;
    301335      timeVariance[i-1] = 0;
    … …  
    305339         
    306340          MCerPhotPix *pix = fEvt->GetPixById(idx);
     341          const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()];
    307342          const MPedestalPix &ped  = (*fPed)[idx];
    308343          const MArrivalTimePix &timepix = (*fTime)[idx];
     344          const Float_t nphot = pix->GetNumPhotons();
    309345
    310346          if (pix == NULL) break;
    … …  
    313349           
    314350            fPixNum[i-1]++;
    315             signal += pix->GetNumPhotons() * (fCam->GetPixRatio(idx));
     351            signal += nphot * (fCam->GetPixRatio(idx));
    316352            noise += pow(ped.GetPedestalRms(),2);
     353           
     354            size += nphot;
     355            meanX += nphot * gpix2.GetX();
     356            meanY += nphot * gpix2.GetY();
    317357           
    318358            time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
    … …  
    337377        }
    338378
    339       // Float_t mean = timeVariance[i-1]/ncore;
    340            
    341       timeVariance[i-1] = 0;
    342      
    343       timeVariance[i-1] = (MAX - MIN)/ncore;
    344       timeVariance[i-1] = (MAX - MIN)/ncore;
     379      meanX /= size;
     380      meanY /= size;
     381 
     382      if (i ==1){
     383        X = meanX;
     384        Y = meanY;
     385      }
     386   
     387      dist = TMath::Power(meanX-X,2) + TMath::Power(meanY-Y,2);
     388      dist = TMath::Sqrt(dist);
     389
     390      timeVariance[i-1] = (MAX - MIN);
    345391
    346392      fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
    … …  
    349395      fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
    350396      fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
    351      
    352     }
     397      fIsl->SetMeanX(i-1, meanX);
     398      fIsl->SetMeanY(i-1, meanY);
     399      fIsl->SetDist(i-1, dist);
     400     
     401    }
    353402 
    354403  fIsl->SetReadyToSave();
    … …  
    493542    }
    494543 
    495  
     544   
    496545  Int_t l = 1;
    497  
     546  Int_t numpixels;
     547  Int_t pixMAX = 0;
     548  Int_t idMAX = 1;
     549
    498550  for(Int_t i = 1;  i<= nvect ; i++)
    499551    {
     552      numpixels = 0;
     553
    500554      if (zeros[i] == 0)
    501555        {
    … …  
    503557            {
    504558              vect[l][k] = vect[i][k];
     559              if (vect[l][k] == 1)
     560                numpixels++;
     561            }
     562          if (numpixels>pixMAX)
     563            {
     564              pixMAX = numpixels;
     565              idMAX = l;
    505566            }
    506567          l++;
    … …  
    509570 
    510571 
     572  //the larger island will correspond to the 1st component of the vector
     573
     574  for(Int_t k = 0; k<npix; k++)
     575    {
     576      vect[nvect+1][k] = vect[1][k];
     577      vect[1][k] = vect[idMAX][k];
     578      vect[idMAX][k] = vect[nvect+1][k];
     579    }
     580 
     581
    511582  //set the number of islands in one event
    512583  fIsl->SetIslNum(fIslNum);
    … …  
    517588  Float_t time[577];
    518589  Float_t timeVariance[fIslNum];
    519  
     590  Float_t size, meanX, meanY, dist;
     591
    520592  //reset the "sets" functions
    521593  if (fIslNum <1)
    … …  
    535607    }
    536608 
     609
     610  Float_t X = 0;
     611  Float_t Y = 0;
    537612   
    538613   for(Int_t i = 1; i<=fIslNum ; i++)
    … …  
    546621      signal = 0;
    547622      noise = 0;
     623
     624      size = 0;
     625      meanX = 0;
     626      meanY = 0;
     627      dist = 0;
     628
    548629      fPixNum[i-1] = 0;
    549630      timeVariance[i-1] = 0;
    … …  
    553634         
    554635          MCerPhotPix *pix = fEvt->GetPixById(idx);
     636          const MGeomPix &gpix2 = (*fCam)[pix->GetPixId()];
    555637          const MPedestalPix &ped  = (*fPed)[idx];
    556638          const MArrivalTimePix &timepix = (*fTime)[idx];
     639          const Float_t nphot = pix->GetNumPhotons();
    557640
    558641          if (pix == NULL) break;
    559642           
    560           if (vect[i][idx] ==1 && fEvt->IsPixelUsed(idx)){
     643          if (vect[i][idx]==1){
    561644           
    562645            fPixNum[i-1]++;
    563             signal += pix->GetNumPhotons() * (fCam->GetPixRatio(idx));
     646            signal += nphot * (fCam->GetPixRatio(idx));
    564647            noise += pow(ped.GetPedestalRms(),2);
     648           
     649            size += nphot;
     650            meanX += nphot * gpix2.GetX();
     651            meanY += nphot * gpix2.GetY();
    565652           
    566653            time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
    … …  
    585672        }
    586673
    587       // Float_t mean = timeVariance[i-1]/ncore;
    588            
    589       timeVariance[i-1] = 0;
    590      
    591       timeVariance[i-1] = (MAX - MIN)/ncore;
    592       timeVariance[i-1] = (MAX - MIN)/ncore;
     674      meanX /= size;
     675      meanY /= size;
     676 
     677      if (i ==1){
     678        X = meanX;
     679        Y = meanY;
     680      }
     681   
     682      dist = TMath::Power(meanX-X,2) + TMath::Power(meanY-Y,2);
     683      dist = TMath::Sqrt(dist);
     684
     685      timeVariance[i-1] = (MAX - MIN);
    593686
    594687      fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
    … …  
    597690      fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
    598691      fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
    599      
    600     }
    601  
     692      fIsl->SetMeanX(i-1, meanX);
     693      fIsl->SetMeanY(i-1, meanY);
     694      fIsl->SetDist(i-1, dist);
     695     
     696    }
     697
     698     
    602699  fIsl->SetReadyToSave();
    603700 

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

    r4328 r4418  
    4747        *fLog << "    + SigToNoise = " << fSigToNoise[i] << endl;
    4848        *fLog << "    + TimeSpread = " << fTimeSpread[i] << endl;
     49        *fLog << "    + MeanX = " << fMeanX[i] << endl;
     50        *fLog << "    + MeanY = " << fMeanY[i] << endl;
     51        *fLog << "    + Dist = " << fDist[i] << endl;
    4952      }
    5053}

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

    r4218 r4418  
    2525    Float_t fTime[20][577];        //  mean of the arrival time 
    2626    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)
    2730
    2831public:
    … …  
    3740    Float_t  GetSigToNoise(Int_t isl)                 { return fSigToNoise[isl]; }
    3841    Float_t  GetArrivalTime(Int_t isl, Int_t idx) { return fTime[isl][idx]; }           
    39     Float_t  GetTimeSpread(Int_t isl) { return fTimeSpread[isl];}         
     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]; }
    4046
    4147    void     SetIslNum(Int_t nisl)                    { fIslNum = nisl; }
    … …  
    4652
    4753    void     SetTimeSpread(Int_t isl, Float_t val) { fTimeSpread[isl] = val; }
    48    
     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; }
     57
     58
    4959    ClassDef(MIslands, 1) // Storage Container for Island Parameters
    5060};
Note: See TracChangeset for help on using the changeset viewer.