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

Timestamp:

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

Author:

aliu

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp/mifae

Files:

• Changelog (modified) (1 diff)
• library/MIslandCalc.cc (modified) (3 diffs)
• library/MIslandCalc.h (modified) (2 diffs)
• programs/makeHillas.cc (modified) (12 diffs)

trunk/MagicSoft/Mars/mtemp/mifae/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2004/05/27 Ester Aliu 
+   * programs/makeHillas.cc
+     - add the possibility of using more than one algorithm to 
+      calculate the islands and use different algorithms for counting
+      islands
+
+   * library/MIslands.[h,cc],MIslandCalc.cc
+     - add a new island cleaning which consists in removing all the
+      islands except the larger one
+     - add a new algorithm of counting islands consisting in consider
+      as the same islands those islands separated by 2 or less pixels
 
  2004/06/02 Javier Rico

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
 }; 

trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc

@@ -65 +65 @@
 
 // initial and final time slices to be used in signal extraction
-const Byte_t hifirst = 0;
-const Byte_t hilast  = 13;
+const Byte_t hifirst = 1;
+const Byte_t hilast  = 14;
 const Byte_t lofirst = 3;
-const Byte_t lolast  = 12;
+const Byte_t lolast  = 14;
 
 // declaration of variables read from datacards
@@ -86 +86 @@
 Float_t  lnew      = 40;
 Int_t    kmethod   = 1;
+Int_t    kalgorithm = 1;
 Int_t    nfiles    = 0;
 
@@ -238 +239 @@
   isl2.SetName("MIslands2");
 
-  if (islflag == 1 || islflag == 2)
+  MIslands      isl3;
+  isl3.SetName("MIslands3");
+
+  if (islflag == 1 || islflag == 2 || islflag == 3)
     {
       plist4.AddToList(&timecam);
@@ -247 +251 @@
     {
       plist4.AddToList(&isl2);
+    }
+
+  if (islflag == 3)
+    {
+      plist4.AddToList(&isl3);
     }
   
@@ -276 +285 @@
   MIslandCalc       island;
   island.SetOutputName("MIslands1");
+  island.SetAlgorithm(kalgorithm);
 
   MBadPixelsTreat   interpolatebadpixels;
@@ -286 +296 @@
   MIslandCalc       island2;
   island2.SetOutputName("MIslands2");  
+  island2.SetAlgorithm(kalgorithm);
+
+  MIslandCalc       island3;
+  island3.SetOutputName("MIslands3");  
   
   
@@ -307 +321 @@
   write->AddContainer("MSrcPosCam"     , "Parameters");
   
-  if (islflag == 1 || islflag == 2)
+  if (islflag == 1 || islflag == 2 || islflag == 3)
     write->AddContainer("MIslands1" , "Parameters");
   if (islflag == 2) 
     write->AddContainer("MIslands2" , "Parameters");
+  if (islflag == 3) 
+    write->AddContainer("MIslands3" , "Parameters");
+
 
   if(display)
@@ -319 +336 @@
       if (islflag == 2)
         disphillas->SetIslandsName("MIslands2");
+      if (islflag == 3)
+        disphillas->SetIslandsName("MIslands3");
     }      
 
@@ -329 +348 @@
   tlist4.AddToList(&clean);
 
-  if (islflag == 1 || islflag == 2)
+  if (islflag == 1 || islflag == 2 || islflag == 3)
     {
       tlist4.AddToList(&timecalc);
@@ -340 +359 @@
       tlist4.AddToList(&island2);
     }
-  
+
+  if (islflag == 3)
+    {
+      tlist4.AddToList(&islclean);
+      tlist4.AddToList(&island3);
+    }
+  
+
   //tlist4.AddToList(&blind2);
   tlist4.AddToList(&hcalc);
@@ -475 +501 @@
         {
           ifun >> islflag;
+
+          // if (islflag == 1 || islflag == 2)
+            ifun >> kalgorithm;
         }
 
@@ -522 +551 @@
   cout << "Cleaning level: ("<<lcore<<","<<ltail<<")" << endl;
   if (islflag == 1 || islflag == 2)
-    cout << "Island calcultation..." << endl;
+    cout << "Island calcultation..." << "using algorithm #" << kalgorithm <<endl;
   if (islflag == 2)
     {