Changeset 4037 for trunk/MagicSoft/Mars/mtemp/MFindStars.cc

Timestamp:

05/10/04 16:53:41 (21 years ago)

Author:

jlopez

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtemp/MFindStars.cc (modified) (16 diffs)

trunk/MagicSoft/Mars/mtemp/MFindStars.cc

@@ -30 +30 @@
 #include "MFindStars.h"
 
+#include <TTimer.h>
+#include <TString.h>
 #include <TFile.h>
 #include <TTree.h>
+#include <TH1F.h>
+#include <TF1.h>
 
 #include "MObservatory.h"
@@ -60 +64 @@
   fTitle = title ? title : "Tool to find stars from DC Currents";
 
+  fNumIntegratedEvents=0;
   fMaxNumIntegratedEvents = 10;
-  fRingInterest = 100.; //[mm] ~ 0.3 deg
-  fMinDCForStars = 3.*fMaxNumIntegratedEvents; //[uA]
+  fRingInterest = 125.; //[mm] ~ 0.4 deg
+  fMinDCForStars = 2.*fMaxNumIntegratedEvents; //[uA]
   
   fPixelsUsed.Set(577);
@@ -166 +171 @@
 Int_t MFindStars::Process()
 {
-    if (fNumIntegratedEvents < fMaxNumIntegratedEvents)
-    {
-      fDisplay.AddCamContent(*fCurr);
-      fNumIntegratedEvents++;
-    }
-    else
+    if (fNumIntegratedEvents >= fMaxNumIntegratedEvents)
     {
       if (fDrive)
@@ -185 +185 @@
       else
         {
+          UInt_t numPixels = fGeomCam->GetNumPixels();
+          TArrayC origPixelsUsed;
+          origPixelsUsed.Set(numPixels);
+
+          for (UInt_t pix=1; pix<numPixels; pix++)
+            {
+              if (fDisplay.IsUsed(pix))
+                origPixelsUsed[pix]=(Char_t)kTRUE;
+              else
+                origPixelsUsed[pix]=(Char_t)kFALSE;
+            }
+          
+          Float_t ped;
+          Float_t rms;
+          DCPedestalCalc(ped, rms);
+          fMinDCForStars = fMinDCForStars>(ped+5*rms)?fMinDCForStars:(ped+5*rms);
+
+          *fLog << dbg << " DC pedestal = " << ped << " pedestal rms = " << rms << endl;
+          *fLog << dbg << " fMinDCForStars " << fMinDCForStars << endl;
+          
+          // Find the star candidats searching the most brights pairs of pixels
           Float_t maxPixelDC;
           MGeomPix maxPixel;
 
-          MHCamera cam = fDisplay;
-          cam.PrintInfo();
-
-          // Find the star candidats searching the most brights pairs of pixels
-         while(1)
+          while(FindPixelWithMaxDC(maxPixelDC, maxPixel))
             {
-              FindPixelWithMaxDC(maxPixelDC, maxPixel);
               *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxPixelDC << " uA) x position(" << setw(3) << maxPixel.GetX() <<  " mm) x position(" << setw(3) << maxPixel.GetY() << " mm)" << endl;
-              
-              if (maxPixelDC<fMinDCForStars)
-                break;
               
               MStarLocalPos *starpos = new MStarLocalPos;
               starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
               starpos->SetCalcValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
+              starpos->SetFitValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
               fStars->GetList()->Add(starpos);
 
@@ -208 +222 @@
               
             }
-         fDisplay = cam;
-         *fLog << dbg <<  GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
+
+          fDisplay.SetUsed(origPixelsUsed);
         }
       
@@ -218 +232 @@
           return kFALSE;
         }
+      else
+          *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
           
       TIter Next(fStars->GetList());
       MStarLocalPos* star;
       while ((star=(MStarLocalPos*)Next()))
-        {
-          FindStar(star);
-          ShadowStar(star);
-        }
-      
+      {
+        FindStar(star);
+        ShadowStar(star);
+      }
       
       //After finding stars reset all vairables
@@ -233 +248 @@
     }
 
+    fDisplay.AddCamContent(*fCurr);
+    fNumIntegratedEvents++;
+
   return kTRUE;
 }
@@ -254 +272 @@
 }
 
+Bool_t MFindStars::DCPedestalCalc(Float_t &ped, Float_t &rms)
+{
+   UInt_t numPixels = fGeomCam->GetNumPixels();
+
+   TH1F dchist("dchist","",120,0.,30.*fMaxNumIntegratedEvents);
+   for (UInt_t pix=1; pix<numPixels; pix++)
+       dchist.Fill(fDisplay.GetBinContent(pix+1));
+
+   Float_t nummaxprobdc = dchist.GetBinContent(dchist.GetMaximumBin());
+   Float_t maxprobdc = dchist.GetBinCenter(dchist.GetMaximumBin());
+   UInt_t bin = dchist.GetMaximumBin();
+   do
+   {
+       bin++;
+   }
+   while(dchist.GetBinContent(bin)/nummaxprobdc > 0.5);
+   Float_t halfmaxprobdc = dchist.GetBinCenter(bin);
+
+   *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
+   *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist.GetBinContent(bin) << "]" << endl;
+
+   Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
+   Float_t min = maxprobdc-3*rmsguess;
+   min = (min<0.?0.:min);
+   Float_t max = maxprobdc+3*rmsguess;
+
+   *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
+
+   TF1 func("func","gaus",min,max);
+   func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
+   
+   dchist.Fit("func","QR0");
+
+   UInt_t aproxnumdegrees = 6*(bin-dchist.GetMaximumBin());
+   Float_t chiq = func.GetChisquare();
+   ped = func.GetParameter(1);
+   rms = func.GetParameter(2);
+
+   *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
+
+   Int_t numsigmas = 5;
+   Axis_t minbin = ped-numsigmas*rms/dchist.GetBinWidth(1);
+   minbin=minbin<1?1:minbin;
+   Axis_t maxbin = ped+numsigmas*rms/dchist.GetBinWidth(1);
+   *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist.Integral((int)minbin,(int)maxbin) << endl;
+   return kTRUE;
+}
+    
 Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix)
 {
@@ -269 +335 @@
             Float_t dc[2];
             dc[0] = fDisplay.GetBinContent(pix+1);
-            
+            if (dc[0] < fMinDCForStars)
+                continue;
+
             MGeomPix &g = (*fGeomCam)[pix];
             Int_t numNextNeighbors = g.GetNumNeighbors();
@@ -280 +348 @@
                     UInt_t swneighbor = g.GetNeighbor(nextNeighbor);
                     dc[1] = fDisplay.GetBinContent(swneighbor+1);
-                    
+                    if (dc[1] < fMinDCForStars)
+                        continue;
                     
                     dcsum = dc[0] + dc[1];
@@ -304 +373 @@
     }
 
+    if (maxDC == 0)
+        return kFALSE;
+
     maxPix = (*fGeomCam)[maxPixIdx[0]];
     return kTRUE;
@@ -311 +383 @@
 {    
 
-  MHCamera cam = fDisplay;
+  UInt_t numPixels = fGeomCam->GetNumPixels();
+  TArrayC origPixelsUsed;
+  origPixelsUsed.Set(numPixels);
   
-    Float_t expX = star->GetXExp();
-    Float_t expY = star->GetYExp();
-
-    UInt_t numPixels = fGeomCam->GetNumPixels();
-
-// First define a area of interest around the expected position of the star
-    for (UInt_t pix=1; pix<numPixels; pix++)
-    {
-
-        Float_t pixXpos=(*fGeomCam)[pix].GetX();
-        Float_t pixYpos=(*fGeomCam)[pix].GetY();
-        Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
-                            (pixYpos-expY)*(pixYpos-expY));
-
-      if (dist < fRingInterest && cam.IsUsed(pix))
-          fPixelsUsed[pix]=(Char_t)kTRUE;
-    }
-
-    cam.SetUsed(fPixelsUsed);
-
-
+  for (UInt_t pix=1; pix<numPixels; pix++)
+    {
+      if (fDisplay.IsUsed(pix))
+        origPixelsUsed[pix]=(Char_t)kTRUE;
+      else
+        origPixelsUsed[pix]=(Char_t)kFALSE;
+    }
+  
+  Float_t expX = star->GetXExp();
+  Float_t expY = star->GetYExp();
+  
+  // First define a area of interest around the expected position of the star
+  for (UInt_t pix=1; pix<numPixels; pix++)
+    {
+      
+      Float_t pixXpos=(*fGeomCam)[pix].GetX();
+      Float_t pixYpos=(*fGeomCam)[pix].GetY();
+      Float_t dist = sqrt((pixXpos-expX)*(pixXpos-expX)+
+                          (pixYpos-expY)*(pixYpos-expY));
+      
+      if ((dist < fRingInterest) && fDisplay.IsUsed(pix))
+        fPixelsUsed[pix]=(Char_t)kTRUE;
+      else
+        fPixelsUsed[pix]=(Char_t)kFALSE;
+    }
+  
+    fDisplay.SetUsed(fPixelsUsed);
+    
     // determine mean x and mean y of the selected px
     Float_t meanX=0;
@@ -343 +424 @@
     for(UInt_t pix=0; pix<numPixels; pix++)
     {
-        if(cam.IsUsed(pix))
+        if(fDisplay.IsUsed(pix))
         {
             usedPx++;
 
-            Float_t charge=cam.GetBinContent(pix+1);
-            Float_t pixXpos=(*fGeomCam)[pix].GetX();
-            Float_t pixYpos=(*fGeomCam)[pix].GetY();
-
-            meanX+=charge*pixXpos;
-            meanY+=charge*pixYpos;
-            meanSqX+=charge*pixXpos*pixXpos;
-            meanSqY+=charge*pixYpos*pixYpos;
-            sumCharge+=charge;
-//          cam.ResetUsed(i); // this px must not be used again!
-        } //for... use  
-    }
-
-    star->SetCalcValues(sumCharge,meanX,meanY,meanSqX,meanSqY);
+            Float_t charge  = fDisplay.GetBinContent(pix+1);
+            Float_t pixXpos = (*fGeomCam)[pix].GetX();
+            Float_t pixYpos = (*fGeomCam)[pix].GetY();
+
+            meanX     += charge*pixXpos;
+            meanY     += charge*pixYpos;
+            meanSqX   += charge*pixXpos*pixXpos;
+            meanSqY   += charge*pixYpos*pixYpos;
+            sumCharge += charge;
+        }
+    }
+
+    *fLog << dbg << " usedPx " << usedPx << endl;
+
+    meanX   /= sumCharge;
+    meanY   /= sumCharge;
+    meanSqX /= sumCharge;
+    meanSqY /= sumCharge;
+
+    Float_t rmsX = TMath::Sqrt(meanSqX - meanX*meanX);
+    Float_t rmsY = TMath::Sqrt(meanSqY - meanY*meanY);
     
+    star->SetCalcValues(sumCharge,meanX,meanY,rmsX,rmsY);
+ 
+    fDisplay.SetUsed(origPixelsUsed);
+
     return kTRUE;
 }
@@ -370 +462 @@
 
 // Define an area around the star which will be set unused.
+    UInt_t shadowPx=0;  
     for (UInt_t pix=1; pix<numPixels; pix++)
     {
-
         Float_t pixXpos  = (*fGeomCam)[pix].GetX();
         Float_t pixYpos  = (*fGeomCam)[pix].GetY();
@@ -378 +470 @@
         Float_t starYpos = star->GetMeanYCalc();
         
-        Float_t starSize = 2*star->GetSigmaMajorAxisCalc();
+        Float_t starSize = 3*star->GetSigmaMajorAxisCalc();
         
         Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
@@ -386 +478 @@
           fPixelsUsed[pix]=(Char_t)kTRUE;
         else
+          {
           fPixelsUsed[pix]=(Char_t)kFALSE;
-    }
+          shadowPx++;
+          }
+    }
+
+    *fLog << dbg << " shadowPx " << shadowPx << endl;
 
     fDisplay.SetUsed(fPixelsUsed);