Changeset 4286 for trunk/MagicSoft/Mars/mtemp/mifae/library

Timestamp:

06/11/04 20:02:08 (21 years ago)

Author:

aliu

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp/mifae/library

Files:

• MIslandCalc.cc (modified) (3 diffs)
• MIslandCalc.h (modified) (2 diffs)

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

@@ -4 +4 @@
 !   Author(s): Ester Aliu, 2/2004 <aliu@ifae.es>
 !   
-!   Last Update: 5/2004
+!   Last Update: 6/2004
 !
 \* ======================================================================== */
@@ -124 +124 @@
 Int_t MIslandCalc::Process()
 {
-
+  
+  if (fIslandAlgorithm == 1)
+    Calc1();
+  
+  if (fIslandAlgorithm == 2)
+    Calc2();
+  return kTRUE;  
+   
+}
+
+
+Int_t MIslandCalc::Calc1(){
+  
+  
+  /////////////////////////////
+  //
+  //        ALGORITHM # 1
+  // counts the number of algorithms as you can see in 
+  // the event display after doing the std image cleaning
+  //
+  /////////////////////////////
+  
   Float_t  noise;
   Float_t  signal;
@@ -333 +354 @@
   fIsl->SetReadyToSave();
   
-  return 1;
-  
+  return kTRUE;  
 }
 
 
+Int_t MIslandCalc::Calc2(){
+  
+  
+  /////////////////////////////
+  //
+  //        ALGORITHM # 2
+  // counts the number of islands considering as the same 
+  // islands the ones separated for 2 or less pixels
+  //
+  /////////////////////////////
+  
+  Float_t  noise;
+  Float_t  signal;
+
+  Int_t    npix = 577;
+
+  Int_t    sflag;
+  Int_t    control;
+  
+  Int_t    nvect = 0;
+  Int_t    fIslNum = 0;
+ 
+  Int_t    vect[50][577];
+
+  Int_t    zeros[50];
+  
+  Int_t kk[577];
+
+  for(Int_t m = 0; m < 50 ; m++)
+    for(Int_t n = 0; n < npix ; n++)
+        vect[m][n] = 0;
+    
+  for(Int_t n = 0; n < 50 ; n++)
+    zeros[n] = 0;
+  
+  for(Int_t n = 0; n < 577 ; n++)
+    kk[n] = 0;
+  
+  MCerPhotPix *pix;
+
+  //1st loop over all pixels
+  MCerPhotEvtIter Next0(fEvt, kFALSE);
+ 
+  while ((pix=static_cast<MCerPhotPix*>(Next0())))
+    {
+      const Int_t idx = pix->GetPixId();
+      
+      const MGeomPix &gpix  = (*fCam)[idx]; 
+      const Int_t    nnmax  = gpix.GetNumNeighbors();
+
+      if( fEvt->IsPixelUsed(idx))
+        {
+          kk[idx] = 1 ;
+          for(Int_t j=0; j< nnmax; j++)
+            {
+              kk[gpix.GetNeighbor(j)] = 1;
+            }
+        } 
+      
+    }
+      
+  //2nd loop over all pixels
+  MCerPhotEvtIter Next(fEvt, kFALSE);
+  
+  while ((pix=static_cast<MCerPhotPix*>(Next())))
+    {
+      const Int_t idx = pix->GetPixId();
+      
+      const MGeomPix &gpix  = (*fCam)[idx];
+      const Int_t    nnmax  = gpix.GetNumNeighbors();
+      
+      if ( kk[idx] > 0)
+        {
+          sflag = 0;
+          
+          for(Int_t j=0; j < nnmax ; j++)
+            {
+              const Int_t idx2 = gpix.GetNeighbor(j);
+              
+              if (idx2 < idx)
+                {
+                  for(Int_t k = 1; k <= nvect; k++)
+                    {
+                      if (vect[k][idx2] == 1)
+                        {
+                          sflag = 1;
+                          vect[k][idx] = 1;
+                        }
+                    }
+                }
+            }
+          
+          if (sflag == 0)
+            {
+              nvect++;
+              vect[nvect][idx] = 1;          
+            }
+          
+        }
+    }
+  
+  fIslNum = nvect;
+  
+  
+  // Repeated Chain Corrections
+  
+  for(Int_t i = 1; i <= nvect; i++)
+    {
+      for(Int_t j = i+1; j <= nvect; j++)
+        {
+          control = 0;
+          for(Int_t k = 0; k < npix; k++)
+            {
+              
+              if (vect[i][k] == 1 && vect[j][k] == 1)
+                {
+                  control = 1; 
+                  break;
+                }
+            }
+          if (control == 1)
+            {
+              for(Int_t k = 0; k < npix; k++)
+                {
+                  if(vect[j][k] == 1)
+                    vect[i][k] = 1;
+                  vect[j][k] = 0;
+                  zeros[j] = 1;
+                }       
+              fIslNum = fIslNum-1;          
+            }
+          
+        }
+    }
+  
+  
+  Int_t l = 1;
+  
+  for(Int_t i = 1;  i<= nvect ; i++)
+    {
+      if (zeros[i] == 0)
+        {
+          for(Int_t k = 0; k<npix; k++)
+            {
+              vect[l][k] = vect[i][k];
+            }
+          l++;
+        }
+    }
+  
+  
+  //set the number of islands in one event
+  fIsl->SetIslNum(fIslNum);
+
+  //examine each island...
+  Int_t fPixNum[fIslNum];
+  Float_t fSigToNoise[fIslNum];
+  Float_t time[577];
+  Float_t timeVariance[fIslNum];
+  
+  //reset the "sets" functions
+  if (fIslNum <1)
+    fIsl->SetIslNum(0);
+
+  for(Int_t i = 0; i<20 ;i++)
+    {
+      fIsl->SetPixNum(i,-1);
+      fIsl->SetSigToNoise(i,-1);
+      fIsl->SetTimeSpread(i,-1);
+    
+      for(Int_t idx = 0; idx<npix; idx++)
+        {
+          fIsl->SetIslId(idx, -1);
+          fIsl->SetArrivalTime(i, idx, -1 );
+        }
+    }
+ 
+   
+   for(Int_t i = 1; i<=fIslNum ; i++)
+    {
+      Int_t n = 0;
+      Int_t ncore = 0;
+      
+      Float_t MIN = 10000;
+      Float_t MAX = 0;
+
+      signal = 0;
+      noise = 0;
+      fPixNum[i-1] = 0;
+      timeVariance[i-1] = 0;
+      
+      for(Int_t idx=0 ; idx<npix ; idx++)
+        {
+          
+          MCerPhotPix *pix = fEvt->GetPixById(idx);
+          const MPedestalPix &ped  = (*fPed)[idx];
+          const MArrivalTimePix &timepix = (*fTime)[idx];
+
+          if (pix == NULL) break;
+            
+          if (vect[i][idx] ==1 && fEvt->IsPixelUsed(idx)){
+            
+            fPixNum[i-1]++;
+            signal += pix->GetNumPhotons() * (fCam->GetPixRatio(idx));
+            noise += pow(ped.GetPedestalRms(),2);
+            
+            time[n] = timepix.IsLoGainUsed() ? timepix.GetArrivalTimeLoGain() : timepix.GetArrivalTimeHiGain();
+            
+            
+            if (fEvt->IsPixelCore(idx)){ 
+              
+              if (time[n] > MAX)
+                MAX = time[n];
+              if (time[n] < MIN)
+                MIN = time[n];
+              
+              ncore++;
+            }
+            
+            fIsl->SetIslId(idx, i-1);
+            fIsl->SetArrivalTime(i-1, idx, time[n]);
+                
+            n++;
+          }
+          
+        }
+
+      // Float_t mean = timeVariance[i-1]/ncore;
+           
+      timeVariance[i-1] = 0;
+     
+      timeVariance[i-1] = (MAX - MIN)/ncore; 
+      timeVariance[i-1] = (MAX - MIN)/ncore; 
+
+      fSigToNoise[i-1]= (Float_t)signal/(Float_t)sqrt(noise);
+
+      fIsl->SetPixNum(i-1,fPixNum[i-1]);
+      fIsl->SetSigToNoise(i-1,fSigToNoise[i-1]);
+      fIsl->SetTimeSpread(i-1, timeVariance[i-1]);
+      
+    } 
+  
+  fIsl->SetReadyToSave();
+  
+  return 1;  
+}

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

@@ -31 +31 @@
   
     TString           fIslName;     //    name of the 'MIslands' container
- 
+  
+    Int_t  fIslandAlgorithm;
+
     Int_t PreProcess(MParList *plist);
     Int_t Process();
+    Int_t Calc1();                  // algorithm of counting islands #1
+    Int_t Calc2();                  // algorithm of counting islands #2    
            
  public:
@@ -39 +43 @@
     void SetOutputName(TString outname)   { fIslName = outname; }
     
+    void SetAlgorithm(Int_t m)   {fIslandAlgorithm = m;}
+    
     ClassDef(MIslandCalc, 0)        // task doing the image cleaning
 };