Changeset 4546

Timestamp:

08/09/04 16:23:06 (20 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mtemp/MFindStars.cc (modified) (33 diffs)
• mtemp/MFindStars.h (modified) (5 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -54 +54 @@
            and MHTelAxisFromStars
        was added
+
+   * mtemp/MFindStars.[h,cc]
+     - add correlated Gauss fit
 
 

trunk/MagicSoft/Mars/mtemp/MFindStars.cc

@@ -15 +15 @@
 ! *
 !
-!   Author(s): Robert Wagner, 2/2004 <mailto:rwagner@mppmu.mpg.de>
-!   Author(s): Javier López , 4/2004 <mailto:jlopez@ifae.es>
+!   Author(s): Robert Wagner,   2/2004 <mailto:rwagner@mppmu.mpg.de>
+!              Javier López ,   4/2004 <mailto:jlopez@ifae.es>
+!              Wolfgang Wittek, 8/2004 <mailto:wittek@mppmu.mpg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2004
@@ -39 +40 @@
 #include <TH1F.h>
 #include <TF1.h>
+#include <TEllipse.h>
+
 
 #include "MObservatory.h"
@@ -70 +73 @@
 //______________________________________________________________________________
 //
-// The 2D gaussian fucntion used to fit the spot of the star
+// The 2D uncorrelated gaussian function used to fit the spot of the star
 //
 static Double_t func(float x,float y,Double_t *par)
 {
-    Double_t value=par[0]*exp(-(x-par[1])*(x-par[1])/(2*par[2]*par[2]))*exp(-(y-par[3])*(y-par[3])/(2*par[4]*par[4]));
+    Double_t value=par[0]*exp( -(x-par[1])*(x-par[1])/(2*par[2]*par[2]) )
+                         *exp( -(y-par[3])*(y-par[3])/(2*par[4]*par[4]) );
     return value;
 }
@@ -80 +84 @@
 //______________________________________________________________________________
 //
+// The 2D correlated gaussian function used to fit the spot of the star
+//
+static Double_t funcCG(float x,float y,Double_t *par)
+{
+  Double_t temp  = 1.0-par[5]*par[5];
+  //  Double_t value = par[0] / (2.0*TMath::Pi()*par[2]*par[4]*sqrt(temp))
+  Double_t value = par[0] 
+    * exp( -0.5/temp * (   (x-par[1])*(x-par[1])/(par[2]*par[2])
+             -2.0*par[5] * (x-par[1])*(y-par[3])/(par[2]*par[4])
+                         + (y-par[3])*(y-par[3])/(par[4]*par[4]) ) );
+    return value;
+}
+
+//______________________________________________________________________________
+//
 // Function used by Minuit to do the fit
 //
@@ -87 +106 @@
   MParList*      plist = (MParList*)gMinuit->GetObjectFit();
   MTaskList*     tlist = (MTaskList*)plist->FindObject("MTaskList");
-  MFindStars*    find = (MFindStars*)tlist->FindObject("MFindStars");
-  MStarLocalCam* stars = (MStarLocalCam*)plist->FindObject("MStarLocalCam");
-  MGeomCam*      geom = (MGeomCam*)plist->FindObject("MGeomCam");
+  MFindStars*    find  = (MFindStars*)tlist->FindObject("MFindStars");
+  MStarLocalCam* stars = 
+          (MStarLocalCam*)plist->FindObject("MStarLocalCam","MStarLocalCam");
+  MGeomCam*      geom  = (MGeomCam*)plist->FindObject("MGeomCam");
 
   MHCamera& display = (MHCamera&)find->GetDisplay();
@@ -100 +120 @@
   UInt_t numPixels = geom->GetNumPixels();
   
+  // fix the correlation if requested
+  if ( find->GetUseCorrelatedGauss() && find->GetFixCorrelation() )
+  {
+    Double_t tandel;
+    if (par[1] != 0.0)
+      tandel = par[3]/par[1];
+    else
+      tandel = 1.e10;
+    Double_t cxx = par[2]*par[2];
+    Double_t cyy = par[4]*par[4];
+    par[5] =   tandel/(tandel*tandel-1.0) * (cyy-cxx)/sqrt(cxx*cyy);
+  }
+
 //calculate chisquare
     Double_t chisq = 0;
@@ -109 +142 @@
     for (UInt_t pixid=1; pixid<numPixels; pixid++) 
     {
+
+
         if (display.IsUsed(pixid))
         {
             x = (*geom)[pixid].GetX();
             y = (*geom)[pixid].GetY();
-            z = display.GetBinContent(pixid+1)-(pixid>lastInnerPixel?outerped:innerped);
+
+            z = display.GetBinContent(pixid+1)
+                -  (pixid>lastInnerPixel?outerped:innerped);
             errorz=(pixid>lastInnerPixel?outerrms:innerrms);
+
 
             if (errorz > 0.0)
             {
               usedPx++;
-              delta  = (z-func(x,y,par))/errorz;
-              chisq += delta*delta;
+
+              Double_t fu;
+              if (find->GetUseCorrelatedGauss())
+                fu = funcCG(x,y,par);
+              else
+                fu = func(x,y,par);
+
+              delta  = (z-fu) / errorz;
+              chisq += delta*delta; 
+
+              if (iflag == 3)
+              {
+                gLog  << "fcn : usedPx, pixid, content, pedestal,z, fu, errorz, delta = "
+                      << usedPx << ",  " << pixid << ",  "
+                      << display.GetBinContent(pixid+1) << ",  "
+                      << (pixid>lastInnerPixel?outerped:innerped)
+                      << ",  " << z << ",  " << fu << ",  " << errorz 
+                      << ",  " << delta << endl;
+              }
             }
             else
@@ -130 +185 @@
     find->SetChisquare(chisq);
     find->SetDegreesofFreedom(usedPx);
+
+    //gLog << "fcn : chisq, usedPx = " << chisq << ",  " << usedPx << endl;
 }
 
+//-------------------------------------------------------------------------
+//
+// Constructor
+//
+
 MFindStars::MFindStars(const char *name, const char *title): 
-  fGeomCam(NULL), fCurr(NULL), fTimeCurr(NULL), fDrive(NULL), fStars(NULL), fNumVar(5)
+  fGeomCam(NULL), fCurr(NULL), fTimeCurr(NULL), fDrive(NULL), fStars(NULL), 
+  fNumVar(6)
 {
   fName  = name  ? name  : "MFindStars";
   fTitle = title ? title : "Tool to find stars from DC Currents";
+
+  // the correlated Gauss function  
+  // is fitted by default
+  fNumVar = 6;
+  fUseCorrelatedGauss = kTRUE;
+  fFixCorrelation     = kFALSE;
 
   fNumIntegratedEvents=0;
@@ -149 +218 @@
   const Float_t pixelSize = 31.5; //[mm]
   fMinuitPrintOutLevel = -1;
+  //fMinuitPrintOutLevel = 3;
   
   fVname = new TString[fNumVar];
@@ -171 +241 @@
   fFix[1]   = 0;
 
-  fVname[2] = "sigmaminor";
+  fVname[2] = "sigmax";
   fVinit[2] = pixelSize;
   fStep[2]  = fVinit[2]/sqrt2;
@@ -185 +255 @@
   fFix[3]   = 0;
 
-  fVname[4] = "sigmamajor";
+  fVname[4] = "sigmay";
   fVinit[4] = pixelSize;
   fStep[4]  = fVinit[4]/sqrt2;
@@ -192 +262 @@
   fFix[4]   = 0;
 
+  if (fUseCorrelatedGauss)
+  {
+    fVname[5] = "xycorr";
+    fVinit[5] = 0.0;
+    fStep[5]  = 0.1;
+    fLimlo[5] = -1.0;
+    fLimup[5] =  1.0;
+    fFix[5]   = 0;
+  }
+
   fObjectFit  = NULL;
   //  fMethod     = "SIMPLEX";
@@ -200 +280 @@
   // Set output level
   //  fLog->SetOutputLevel(3); // No dbg messages
+
+  fGeometryFile="";
+  fBSCFile="";
 }
+
+//-------------------------------------------------------------------------
+//
+// Set 'fUseCorrelatedGauss' and 'fFixCorrelation'  flag
+//
+//     if 'fUseCorrelatedGauss' is TRUE a 2dim Gaussian with correlation
+//                                 will be fitted
+//     if 'fFixCorrelation' is TRUE the orientation of the 2dim Gaussian
+//                             will be kept fixed at the angle delta,
+//                             where tan(delta) = ymean/xmean
+
+void MFindStars::SetUseCorrelatedGauss(Bool_t usecorrgauss, Bool_t fixcorr)
+{
+  fUseCorrelatedGauss = usecorrgauss;
+  fFixCorrelation   = fixcorr;
+
+  if (usecorrgauss)
+    fNumVar = 6;
+  else
+    fNumVar = 5;
+}
+
+//-------------------------------------------------------------------------
+//
+// PreProcess
+//
 
 Int_t MFindStars::PreProcess(MParList *pList)
@@ -233 +342 @@
     }
 
-    //    fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
+    fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
 
     if (!fDrive)
       {
+
         *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl;
+
       }
     else
       {
-//         //Initialitation MAstroCamera
-//         // Name of a MC file having MGeomCam and MMcConfigRunHeader
-//         TString fname = "../data/Gamma_zbin9_90_7_1480to1489_w0.root";
-
-//         // Time for which to get the picture
-//         //    MTime time;
-//         //    time.Set(2004, 2, 28, 01, 32, 15);
-        
-//         // Current observatory
-//         MObservatory magic1;
-        
-//         // Right Ascension [h] and declination [deg] of source
-//         // Currently 'perfect' pointing is assumed
-//         //    const Double_t ra  = MAstro::Hms2Rad(5, 34, 31.9);
-//         //    const Double_t dec = MAstro::Dms2Rad(22, 0, 52.0);
-
-//         // --------------------------------------------------------------------------
-//         // Create camera display from geometry
-//         MMcConfigRunHeader *config=0;
-//         MGeomCam           *geom=0;
-        
-//         TFile file(fname);
-//         TTree *tree = (TTree*)file.Get("RunHeaders");
-//         tree->SetBranchAddress("MMcConfigRunHeader", &config);
-//         if (tree->GetBranch("MGeomCam"))
-//           tree->SetBranchAddress("MGeomCam", &geom);
-//         tree->GetEntry(0);
-        
-//         fAstro.SetMirrors(*config->GetMirrors());
-//         fAstro.SetGeom(*geom);
-        
-        
-//         fAstro.SetLimMag(6);
-//         fAstro.SetRadiusFOV(3);
-//         fAstro.ReadBSC("../data/bsc5.dat");
-        
-//         fAstro.SetObservatory(magic1);
-//         // Time for which to get the picture
-//         //    MTime time;
-//         //    time.Set(2004, 2, 28, 01, 32, 15);
-//         //   fAstro.SetTime(time);
-//         fAstro.SetTime(*fTimeCurr);
-//         fAstro.SetGuiActive();
-        
-//         fAstro.SetStarList(fStars->GetList());
-        
+        
+        MObservatory magic1;
+
+        MMcConfigRunHeader *config=0;
+        MGeomCam           *geom=0;
+
+        TFile file(fGeometryFile);
+        TTree *tree = (TTree*)file.Get("RunHeaders");
+        tree->SetBranchAddress("MMcConfigRunHeader", &config);
+        if (tree->GetBranch("MGeomCam")) tree->SetBranchAddress("MGeomCam", &geom);
+        tree->GetEntry(0);
+        
+        fAstro.SetMirrors(*config->GetMirrors());
+        fAstro.SetGeom(*geom);  
+        fAstro.ReadBSC(fBSCFile);
+        
+        fAstro.SetObservatory(magic1);
+        
       }
     
 
-    fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"));
+    fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam"),"MStarLocalCam");
     if (!fStars)
     {
@@ -312 +395 @@
     gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
 
+    // suppress warnings
+    Int_t errWarn;
+    Double_t tmpwar = 0;
+    gMinuit->mnexcm("SET NOW", &tmpwar, 0, errWarn);
+
+
     // Set MParList object pointer to allow mimuit to access internally
     gMinuit->SetObjectFit(pList);
@@ -317 +406 @@
     return kTRUE;
 }
+
+//------------------------------------------------------------------------
+//
+// Process :
+//              
+//
+//
+//
 
 Int_t MFindStars::Process()
@@ -325 +422 @@
   origPixelsUsed.Set(numPixels);
 
-    if (fNumIntegratedEvents >= fMaxNumIntegratedEvents)
-    {
-
-      if (fDrive)
-        {
-//           Float_t ra  = fDrive->GetRa();
-//           Float_t dec = fDrive->GetDec();
+  if (fNumIntegratedEvents >= fMaxNumIntegratedEvents) {
+
+    //Fist delete the previous stars in the list
+    fStars->GetList()->Delete();
+
+    if (fDrive) {
+
+      //If drive information is provided we take RaDec info
+      //from the drive and let the star list fill by the astrocamera.
+
+      //FIXME: rwagner: Doesn't work as expected
+      //FIXME: rwagner: For the time being set manually.
+      //fAstro.SetRaDec(fDrive->GetRa(), fDrive->GetDec());              
+
+      fAstro.SetTime(*fTimeCurr);
+      fAstro.SetGuiActive();
+      fAstro.FillStarList(fStars->GetList());      
+
+      cout << "Number of Stars added by astrocamera is " <<fStars->GetList()->GetSize() << endl;
+      
+      MStarLocalPos* starpos;
+      TIter Next(fStars->GetList());
+      while ((starpos=(MStarLocalPos*)Next())) {
+        starpos->SetCalcValues(40,40,starpos->GetXExp(),starpos->GetYExp(),fRingInterest/2,fRingInterest/2);
+        starpos->SetFitValues (40,40,starpos->GetXExp(),starpos->GetYExp(),fRingInterest/2,fRingInterest/2,0.,-1);
+      }
+
+      for (UInt_t pix=1; pix<numPixels; pix++) {
+        if (fDisplay.IsUsed(pix))
+          origPixelsUsed[pix]=(Char_t)kTRUE;
+        else
+          origPixelsUsed[pix]=(Char_t)kFALSE;
+      }
+      
+      DCPedestalCalc();
+
+    } 
+    else 
+    {
+      
+      cout << "No drive information available: Will not use a star catalog to identify stars."<< endl;
+      
+      for (UInt_t pix=1; pix<numPixels; pix++) {
+        if (fDisplay.IsUsed(pix))
+          origPixelsUsed[pix]=(Char_t)kTRUE;
+        else
+          origPixelsUsed[pix]=(Char_t)kFALSE;
+      }
           
-//           fAstro.SetRaDec(ra, dec);
-//           fAstro.SetGuiActive();
-          
-//           fAstro.FillStarList();
-        }
+      if (DCPedestalCalc()) {
+
+        Float_t innermin = fStars->GetInnerPedestalDC()+fDCTailCut*fStars->GetInnerPedestalRMSDC();
+        Float_t outermin = fStars->GetOuterPedestalDC()+fDCTailCut*fStars->GetOuterPedestalRMSDC();
+        fMinDCForStars = innermin>outermin?innermin:outermin;
+        if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars = " << fMinDCForStars << endl;
+              
+        // Find the star candidates searching the most brights pairs of pixels
+        Float_t maxPixelDC;
+        MGeomPix maxPixel;
+        
+        while(FindPixelWithMaxDC(maxPixelDC, maxPixel)) {
+          
+          MStarLocalPos *starpos = new MStarLocalPos;
+          starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
+          starpos->SetCalcValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
+          starpos->SetFitValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2,0.,-1);
+          fStars->GetList()->Add(starpos);
+          
+          ShadowStar(starpos);
+        }       
+        fDisplay.SetUsed(origPixelsUsed);
+      }      
+    }
+
+    //Show the stars found
+    //fStars->Print("namepossizechierr");
+   
+    //loop to extract position of stars on the camera
+    if (fStars->GetList()->GetSize() == 0) {
+      *fLog << err << GetName() << "No stars candidates in the camera." << endl;
+      return kCONTINUE;
+    } else
+      *fLog << inf << GetName() << " found " << fStars->GetList()->GetSize() 
+            << " star candidates in the camera." << endl;
+    
+    for (UInt_t pix=1; pix<numPixels; pix++) {
+      if (fDisplay.IsUsed(pix))
+        origPixelsUsed[pix]=(Char_t)kTRUE;
       else
-        {
-          //Fist delete the previus stars in the list
-          fStars->GetList()->Delete();
-          //
-
-          for (UInt_t pix=1; pix<numPixels; pix++)
-            {
-              if (fDisplay.IsUsed(pix))
-                origPixelsUsed[pix]=(Char_t)kTRUE;
-              else
-                  origPixelsUsed[pix]=(Char_t)kFALSE;
-            }
-          
-          if (DCPedestalCalc())
-            {
-
-              Float_t innermin = fStars->GetInnerPedestalDC()+fDCTailCut*fStars->GetInnerPedestalRMSDC();
-              Float_t outermin = fStars->GetOuterPedestalDC()+fDCTailCut*fStars->GetOuterPedestalRMSDC();
-              fMinDCForStars = innermin>outermin?innermin:outermin;
-              if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars = " << fMinDCForStars << endl;
-              
-              // Find the star candidats searching the most brights pairs of pixels
-              Float_t maxPixelDC;
-              MGeomPix maxPixel;
-
-              while(FindPixelWithMaxDC(maxPixelDC, maxPixel))
-                {
-                  
-                  MStarLocalPos *starpos = new MStarLocalPos;
-                  starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY());
-                  starpos->SetCalcValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2);
-                  starpos->SetFitValues(maxPixelDC,maxPixelDC,maxPixel.GetX(),maxPixel.GetY(),fRingInterest/2,fRingInterest/2,0.,1);
-                  fStars->GetList()->Add(starpos);
-                  
-                  ShadowStar(starpos);
-                  
-                }
-              
-              fDisplay.SetUsed(origPixelsUsed);
-            }
-          
-        }
-      
-      //loop to extract position of stars on the camera
-      if (fStars->GetList()->GetSize() == 0)
-        {
-          *fLog << err << GetName() << " No stars candidates in the camera." << endl;
-          return kCONTINUE;
-        }
-      else
-          *fLog << inf << GetName() << " Found " << fStars->GetList()->GetSize() << " stars candidates in the camera." << endl;
-          
-
-      for (UInt_t pix=1; pix<numPixels; pix++)
-        {
-          if (fDisplay.IsUsed(pix))
-            origPixelsUsed[pix]=(Char_t)kTRUE;
-          else
-            origPixelsUsed[pix]=(Char_t)kFALSE;
-        }
-
-      TIter Next(fStars->GetList());
-      MStarLocalPos* star;
-      while ((star=(MStarLocalPos*)Next()))
-        {
-          FindStar(star);
-          ShadowStar(star);
-        }
-      
-      //After finding stars reset all vairables
-      fDisplay.Reset();
-      fStars->GetDisplay().Reset(); //FIXME: Put this display just in the container
-      fDisplay.SetUsed(origPixelsUsed);
-      fNumIntegratedEvents=0;
-    }
-
-    for (UInt_t pix=1; pix<numPixels; pix++)
-      {
-        if (fDisplay.IsUsed(pix))
-          origPixelsUsed[pix]=(Char_t)kTRUE;
-        else
-          origPixelsUsed[pix]=(Char_t)kFALSE;
-        
-      }
-
-    fDisplay.AddCamContent(*fCurr);
-    fNumIntegratedEvents++;
+        origPixelsUsed[pix]=(Char_t)kFALSE;
+    }
+
+    TIter Next(fStars->GetList());
+    MStarLocalPos* star;
+    while ((star=(MStarLocalPos*)Next())) {
+       FindStar(star);
+       ShadowStar(star);
+    }
+
+    //Show the stars found: Here it is interesting to see what FindStars
+    //made out of the positions we gave to it.
+    if (fMinuitPrintOutLevel>=0)
+      fStars->Print("namepossizchierr");
+
+    //After finding stars reset all variables
+    fDisplay.Reset();
+    fStars->GetDisplay().Reset(); //FIXME: Put this display just in the container
     fDisplay.SetUsed(origPixelsUsed);
+    fNumIntegratedEvents=0;
+  }
+
+  for (UInt_t pix=1; pix<numPixels; pix++) {
+    if (fDisplay.IsUsed(pix))
+      origPixelsUsed[pix]=(Char_t)kTRUE;
+    else
+      origPixelsUsed[pix]=(Char_t)kFALSE;
     
-
+  }
+  
+  fDisplay.AddCamContent(*fCurr);
+  fNumIntegratedEvents++;
+  fDisplay.SetUsed(origPixelsUsed);
+  
   return kTRUE;
 }
@@ -449 +565 @@
     fDisplay.SetUsed(fPixelsUsed);
 }
+
+//-------------------------------------------------------------------------
+//
+// DCPedestalCalc :
+//
+//  
+//
+//
 
 Bool_t MFindStars::DCPedestalCalc()
@@ -463 +587 @@
    Float_t rms;
 
-   TH1F **dchist = new TH1F*[2];
+   //TH1F **dchist = new TH1F*[2];
+   TH1F *dchist[2];
+   TString tit;
+   TString nam;
    for (UInt_t i=0; i<2; i++)
-      dchist[i] = new TH1F("","",26,0.4*fMaxNumIntegratedEvents,3.*fMaxNumIntegratedEvents);
-   
+   {
+      nam = i;
+      tit = "dchist[";
+      tit += i;
+      tit += "]";
+      dchist[i] = new TH1F(nam, tit,
+                  26,0.4*fMaxNumIntegratedEvents,3.*fMaxNumIntegratedEvents);
+   }   
+
    for (UInt_t pix=1; pix<=lastInnerPixel; pix++)
        dchist[0]->Fill(fDisplay.GetBinContent(pix+1));
@@ -474 +608 @@
    // inner/outer pixels
    for (UInt_t i=0; i<2; i++)
-    {
+   {
       Float_t nummaxprobdc = dchist[i]->GetBinContent(dchist[i]->GetMaximumBin());
       Float_t maxprobdc = dchist[i]->GetBinCenter(dchist[i]->GetMaximumBin());
@@ -517 +651 @@
         {
           *fLog << dbg << "dchist[i]->GetEntries()" << dchist[i]->GetEntries();
+//            TCanvas *c1 = new TCanvas("c2","c2",500,800);
+//            dchist[i]->Draw();
+//            c1->Print("dchist.ps");
+//            delete c1;
+//            exit(1);
         }
       else if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
@@ -527 +666 @@
         }
    
-//       TCanvas *c1 = new TCanvas("c2","c2",500,800);
-//       dchist[i]->Draw();
-//       c1->Print("dchist.ps");
-//       getchar();
-//       delete c1;
-//       exit(1);
-
       if (i == 0)
         {
@@ -544 +676 @@
           fStars->SetOuterPedestalRMSDC(rms);
         }
-
-      
-
-    }
+   }
    
 
+
    for (UInt_t i=0; i<2; i++)
+   {
       delete dchist[i];
-   delete [] dchist;
+   }
+   //delete [] dchist;
 
    //=================================================================
@@ -632 +764 @@
 }
 
+
+//----------------------------------------------------------------------------
+//
+// FindStar
+//
+//
+//
+
 Bool_t MFindStars::FindStar(MStarLocalPos* star)
 {    
-
   UInt_t numPixels = fGeomCam->GetNumPixels();
   Float_t innerped = fStars->GetInnerPedestalDC();
@@ -652 +791 @@
   Float_t expX = star->GetXExp();
   Float_t expY = star->GetYExp();
-  
+
   Float_t max=0;
   UInt_t  pixmax=0;
@@ -660 +799 @@
   Float_t meanSqY=0;
   Float_t sumCharge=0;
-  UInt_t usedInnerPx=0; 
-  UInt_t usedOuterPx=0; 
+  UInt_t  usedInnerPx=0;        
+  UInt_t  usedOuterPx=0;        
+
+  Float_t meanXold = 1.e10;
+  Float_t meanYold = 1.e10;
 
   const Float_t meanPresition = 3.; //[mm]
@@ -667 +809 @@
   UInt_t numIterations = 0;
 
-  do
-    {
-  // 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);
+  //--------------------   start of iteration loop   -----------------------
+  for (UInt_t it=0; it<maxNumIterations; it++)
+  {
+      //rwagner: Need to find px with highest dc and need to assume it is
+      // somewhere near the center of the star
+      //after that we need to REDEFINE our roi! because expected pos could be
+      // quite off that px!
+      
+      // 1.) Initial roi around expected position
+
+      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);
+
+      // 2.) Find px with highest dc in that region
+
+      max    = 0.0;
+      pixmax =   0;
+      for(UInt_t pix=0; pix<numPixels; pix++)   
+        if(fDisplay.IsUsed(pix))
+          {
+            Float_t charge  = fDisplay.GetBinContent(pix+1);          
+            if (charge>max)
+              {
+                max=charge;
+                pixmax=pix;
+              }
+          }         
+
+      // 3.) make it new center
+
+      expX = (*fGeomCam)[pixmax].GetX();
+      expY = (*fGeomCam)[pixmax].GetY();
+
+      for (UInt_t pix=1; pix<numPixels; pix++)
+        fPixelsUsed[pix]=(Char_t)kTRUE;       
+      fDisplay.SetUsed(fPixelsUsed);
+
+      // 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
-    usedInnerPx=0;      
-    usedOuterPx=0;      
-    for(UInt_t pix=0; pix<numPixels; pix++)
-    {
-        if(fDisplay.IsUsed(pix))
+      // determine mean x and mean y
+      max         = 0.0;
+      pixmax      =   0;
+      meanX       = 0.0;
+      meanY       = 0.0;
+      meanSqX     = 0.0;
+      meanSqY     = 0.0;
+      sumCharge   = 0.0;
+      usedInnerPx =   0;        
+      usedOuterPx =   0;        
+
+      for(UInt_t pix=0; pix<numPixels; pix++)
         {
-            pix>lastInnerPixel?usedOuterPx++:usedInnerPx++;
-
-            Float_t charge  = fDisplay.GetBinContent(pix+1);
-            Float_t pixXpos = (*fGeomCam)[pix].GetX();
-            Float_t pixYpos = (*fGeomCam)[pix].GetY();
-
-            if (charge>max)
-              {
-                max=charge;
-                pixmax=pix;
-              }
-            
-            meanX     += charge*pixXpos;
-            meanY     += charge*pixYpos;
-            meanSqX   += charge*pixXpos*pixXpos;
-            meanSqY   += charge*pixYpos*pixYpos;
-            sumCharge += charge;
-        }
-    }
-
-    if (fMinuitPrintOutLevel>=0) *fLog << dbg << " usedInnerPx " << usedInnerPx << " usedOuterPx " << usedOuterPx << endl;
-
-    meanX   /= sumCharge;
-    meanY   /= sumCharge;
-    meanSqX /= sumCharge;
-    meanSqY /= sumCharge;
-
-    expX = meanX;
-    expY = meanY;
-
-    if (++numIterations >  maxNumIterations)
-      {
-        *fLog << warn << GetName() << "Mean calculation not converge after " << maxNumIterations << " iterations" << endl;
+          if(fDisplay.IsUsed(pix))
+            {
+              pix>lastInnerPixel?usedOuterPx++:usedInnerPx++;
+              
+              Float_t charge  = fDisplay.GetBinContent(pix+1);
+              Float_t pixXpos = (*fGeomCam)[pix].GetX();
+              Float_t pixYpos = (*fGeomCam)[pix].GetY();
+              
+              if (charge>max)
+                {
+                  max=charge;
+                  pixmax=pix;
+                }
+              
+              meanX     += charge*pixXpos;
+              meanY     += charge*pixYpos;
+              meanSqX   += charge*pixXpos*pixXpos;
+              meanSqY   += charge*pixYpos*pixYpos;
+              sumCharge += charge;
+            }
+        } 
+      
+      if (fMinuitPrintOutLevel>=0) *fLog << dbg << " usedInnerPx " << usedInnerPx << " usedOuterPx " << usedOuterPx << endl;
+      
+      meanX   /= sumCharge;
+      meanY   /= sumCharge;
+      meanSqX /= sumCharge;
+      meanSqY /= sumCharge;
+
+      // stop iteration when (meanX, meanY) becomes stable
+      if ( TMath::Abs(meanX-meanXold) < meanPresition  && 
+           TMath::Abs(meanY-meanYold) < meanPresition    )
         break;
-      }
-        
-    }while(TMath::Abs(meanX-expX) > meanPresition || TMath::Abs(meanY-expY) > meanPresition);
-  
-    Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12;
-    Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12;
-
-    if ( rmsX > 0)
-      rmsX =  TMath::Sqrt(rmsX);
-    else
-      {
-        *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl;
-        *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl;
-        rmsX = 0.;
-      }
-
-    if ( rmsY > 0)
-      rmsY =  TMath::Sqrt(rmsY);
-    else
-      {
-        *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl;
-        *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl;
-        rmsY = 0.;
-      }
-    
-    // Substrack pedestal DC
-    sumCharge-= (usedInnerPx*innerped+usedOuterPx*outerped)/(usedInnerPx+usedOuterPx);
-    max-=pixmax>lastInnerPixel?outerped:innerped;
-    
-
-    star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY);
-
-    if (rmsX <= 0. || rmsY <= 0.)
-      return kFALSE;
+
+      meanXold = meanX;
+      meanYold = meanY;
+      numIterations++;
+  }
+  // this statement was commented out because the condition 
+  // was never fulfilled :
+  //while(TMath::Abs(meanX-expX) > meanPresition || 
+  //      TMath::Abs(meanY-expY) > meanPresition);
+  //--------------------   end of iteration loop   -----------------------
+
+  if (numIterations >=  maxNumIterations)
+  {
+     *fLog << warn << GetName() << "Mean calculation not converged after " 
+           << maxNumIterations << " iterations" << endl;
+  }
+  
+
+  //Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12;
+  //Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12;
+
+  Float_t rmsX = (meanSqX - meanX*meanX);
+  Float_t rmsY = (meanSqY - meanY*meanY);
+  
+  if ( rmsX > 0)
+    rmsX =  TMath::Sqrt(rmsX);
+  else
+    {
+      *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl;
+      *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl;
+      rmsX = 0.;
+    }
+  
+  if ( rmsY > 0)
+    rmsY =  TMath::Sqrt(rmsY);
+  else
+    {
+      *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl;
+      *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl;
+      rmsY = 0.;
+    }
+  
+  // Substrack pedestal DC
+  sumCharge-= (usedInnerPx*innerped+usedOuterPx*outerped)/(usedInnerPx+usedOuterPx);
+  max-=pixmax>lastInnerPixel?outerped:innerped;
+  
+  
+  star->SetCalcValues( sumCharge, max, meanX, meanY, rmsX, rmsY);
+  
+  if (rmsX <= 0. || rmsY <= 0.)
+    return kFALSE;
     
     
@@ -772 +982 @@
     if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars " << fMinDCForStars << " pixmax>lastInnerPixel?outerped:innerped " << (pixmax>lastInnerPixel?outerped:innerped) << " fMaxNumIntegratedEvents " << fMaxNumIntegratedEvents << endl;
 
-  //Initialate variables for fit
+  //Initialize variables for fit
     fVinit[0] = max;
     fLimlo[0] = fMinDCForStars-(pixmax>lastInnerPixel?outerped:innerped);
@@ -782 +992 @@
     //Init steps
     for(Int_t i=0; i<fNumVar; i++)
-      {
+    {
         if (fVinit[i] != 0)
           fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
+        else
+          fStep[i] = 0.1;
+
         if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fVinit[" << i << "] " << fVinit[i];
         if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fStep[" << i << "] " << fStep[i];
         if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimlo[" << i << "] " << fLimlo[i];
         if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimup[" << i << "] " << fLimup[i] << endl;
-      }
-    //
+    }
+
+    // set the correlation fixed if requested
+    if (fUseCorrelatedGauss && fFixCorrelation)
+    {
+      fStep[5] = 0.0;
+      fFix[5]  = 1;
+    }
+
 
     // -------------------------------------------
@@ -809 +1029 @@
     gMinuit->mnexcm(fMethod, arglist ,2,ierflg);
 
+    // call fcn with iflag = 3 
+    //arglist[0] = 3;
+    //Int_t ierflg3;
+    //gMinuit->mnexcm("CALL", arglist, 1, ierflg3);     
+
     clock.Stop();
 
@@ -817 +1042 @@
       }
 
-    Double_t integratedCharge;
-    Double_t maxFit, maxFitError;
-    Double_t meanXFit, meanXFitError;
-    Double_t sigmaMinorAxis, sigmaMinorAxisError;
-    Double_t meanYFit, meanYFitError;
-    Double_t sigmaMajorAxis, sigmaMajorAxisError;
-    Float_t chisquare = GetChisquare();
-    Int_t   dregrees  = GetDegreesofFreedom()-fNumVar;
-
-    if (!ierflg)
+
+    //----------  for the uncorrelated Gauss fit (start)   ---------
+    if (!fUseCorrelatedGauss)
+    {
+      Double_t integratedCharge;
+      Double_t maxFit,   maxFitError;
+      Double_t meanXFit, meanXFitError;
+      Double_t sigmaX,   sigmaXError;
+      Double_t meanYFit, meanYFitError;
+      Double_t sigmaY,   sigmaYError;
+      Float_t  chisquare = GetChisquare();
+      Int_t    dregrees  = GetDegreesofFreedom()-fNumVar;
+
+      if (!ierflg)
       {
-        gMinuit->GetParameter(0,maxFit, maxFitError);
-        gMinuit->GetParameter(1,meanXFit,meanXFitError);
-        gMinuit->GetParameter(2,sigmaMinorAxis,sigmaMinorAxisError);
-        gMinuit->GetParameter(3,meanYFit,meanYFitError);
-        gMinuit->GetParameter(4,sigmaMajorAxis,sigmaMajorAxisError);
+        gMinuit->GetParameter(0, maxFit,   maxFitError);
+        gMinuit->GetParameter(1, meanXFit, meanXFitError);
+        gMinuit->GetParameter(2, sigmaX,   sigmaXError);
+        gMinuit->GetParameter(3, meanYFit, meanYFitError);
+        gMinuit->GetParameter(4, sigmaY,   sigmaYError);
         
         //FIXME: Do the integral properlly
         integratedCharge = 0.;
-
-        
-      }
-    else
+      }
+      else
       {
         maxFit = 0.;
         meanXFit = 0.;
-        sigmaMinorAxis = 0.;
+        sigmaX = 0.;
         meanYFit = 0.;
-        sigmaMajorAxis = 0.;
+        sigmaY = 0.;
         integratedCharge = 0.;
 
         *fLog << err << "TMinuit::Call error " << ierflg << endl;
       }
-    
-    star->SetFitValues(integratedCharge,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,chisquare,dregrees);
-    
+
+      //rwagner: get error matrix
+      Double_t matrix[5][5];
+      gMinuit->mnemat(&matrix[0][0],5);
+
+      star->SetFitValues(integratedCharge,maxFit,       meanXFit,    meanYFit,
+                         sigmaX,          sigmaY,       chisquare,   dregrees,
+                         matrix[1][1],    matrix[1][3], matrix[3][3]);
+
+      // set the results from the correlated Gauss fit to zero
+      star->SetCGFitValues(100.0,  100.0,  0.0, 0.0,  0.0,   0.0,  0.0,
+                           0.0,      0.0,  0.0, 0.0,   -1);
+    }
+    //----------  for the uncorrelated Gauss fit (end)   ---------    
+
+    //----------  for the correlated Gauss fit (start)   ---------
+    if (fUseCorrelatedGauss)
+    {
+      TArrayD fValues(fNumVar);
+      TArrayD fErrors(fNumVar);
+
+      const static Int_t fNdim = 6;
+      Double_t matrix[fNdim][fNdim];
+
+      Double_t fmin   = 0.0;
+      Double_t fedm   = 0.0;
+      Double_t errdef = 0.0;
+      Int_t npari     =   0;
+      Int_t nparx     =   0;
+      Int_t istat     =   0;
+      gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
+      if (fMinuitPrintOutLevel>=0)
+      { 
+            *fLog << dbg  
+            << "Status of Correlated Gauss fit to the DC currents : " << endl;
+      }
+      if (fMinuitPrintOutLevel>=0) 
+      {
+            *fLog << dbg 
+            << "       fmin, fedm, errdef, npari, nparx, istat = "
+            << fmin << ",  " << fedm << ",  " << errdef << ",  " << npari
+            << ",  " << nparx << ",  " << istat << endl;
+      }
+
+      if (!ierflg)
+      {
+        for (Int_t k=0; k<fNumVar; k++)
+          gMinuit->GetParameter( k, fValues[k], fErrors[k] );
+      }
+      else
+      {
+        *fLog << err << "Correlated Gauss fit failed : ierflg, istat = " 
+              << ierflg << ",  " << istat << endl;
+      }
+
+      Float_t  charge = 0.0;
+      Float_t  chi2   = fmin;
+      Int_t    ndof   = GetDegreesofFreedom()-fNumVar;
+
+      //get error matrix
+      if (istat >= 1)
+      {
+        gMinuit->mnemat( &matrix[0][0], fNdim);
+
+        // copy covariance matrix into a matrix which includes also the fixed
+        // parameters
+        TString  name;
+        Double_t bnd1, bnd2, val, err;
+        Int_t    jvarbl;
+        Int_t    kvarbl;
+        for (Int_t j=0; j<fNumVar; j++)
+        {
+            if (gMinuit)
+                gMinuit->mnpout(j, name, val, err, bnd1, bnd2, jvarbl);
+
+            for (Int_t k=0; k<fNumVar; k++)
+            {
+              if (gMinuit)
+                  gMinuit->mnpout(k, name, val, err, bnd1, bnd2, kvarbl);
+
+              matrix[j][k] = jvarbl==0 || kvarbl==0 ? 0 : matrix[jvarbl-1][kvarbl-1];
+            }
+        }
+      }
+      else
+      {
+        // error matrix was not calculated, construct it
+        if (fMinuitPrintOutLevel>=0)
+        {
+        *fLog << warn 
+              << "error matrix is not defined, construct it" << endl;
+        }
+        for (Int_t j=0; j<fNumVar; j++)
+        {
+          for (Int_t k=0; k<fNumVar; k++)
+              matrix[j][k] = 0;
+
+          matrix[j][j] = fErrors[j]*fErrors[j];
+        }
+      }
+
+      if (fMinuitPrintOutLevel>=0)
+      {
+        *fLog << "Before calling SetCGFitValues : " << endl; 
+        *fLog << "fValues = " << fValues[0] << ",  " << fValues[1] << ",  "
+              << fValues[2] << ",  " << fValues[3] << ",  " << fValues[4]
+              << ",  " << fValues[5] << endl;
+
+        *fLog << "fErrors = " << fErrors[0] << ",  " << fErrors[1] << ",  "
+              << fErrors[2] << ",  " << fErrors[3] << ",  " << fErrors[4]
+              << ",  " << fErrors[5] << endl;
+
+        *fLog << "error matrix = " << matrix[0][0] << "   " << matrix[0][1]
+                          << "   " << matrix[0][2] << "   " << matrix[0][3]
+                          << "   " << matrix[0][4] << "   " << matrix[0][5]
+              << endl;
+        *fLog << "               " << matrix[1][0] << "   " << matrix[1][1]
+                          << "   " << matrix[1][2] << "   " << matrix[1][3]
+                          << "   " << matrix[1][4] << "   " << matrix[1][5]
+              << endl;
+        *fLog << "               " << matrix[2][0] << "   " << matrix[2][1]
+                          << "   " << matrix[2][2] << "   " << matrix[2][3]
+                          << "   " << matrix[2][4] << "   " << matrix[2][5]
+              << endl;
+        *fLog << "               " << matrix[3][0] << "   " << matrix[3][1]
+                          << "   " << matrix[3][2] << "   " << matrix[3][3]
+                          << "   " << matrix[3][4] << "   " << matrix[3][5]
+              << endl;
+        *fLog << "               " << matrix[4][0] << "   " << matrix[4][1]
+                          << "   " << matrix[4][2] << "   " << matrix[4][3]
+                          << "   " << matrix[4][4] << "   " << matrix[4][5]
+              << endl;
+        *fLog << "               " << matrix[5][0] << "   " << matrix[5][1]
+                          << "   " << matrix[5][2] << "   " << matrix[5][3]
+                          << "   " << matrix[5][4] << "   " << matrix[5][5]
+              << endl;
+      }
+
+      star->SetCGFitValues(charge,       fValues[0],  fValues[1], fValues[3],
+                           fValues[2],   fValues[4],  fValues[5],
+                           matrix[1][1], matrix[1][3],matrix[3][3],
+                           chi2,         ndof);
+      // set the results from the uncorrelated Gauss fit to zero
+      star->SetFitValues(100.0, 100.0,  0.0, 0.0,  0.0, 0.0,  0.0, -1,
+                           0.0, 0.0, 0.0);
+    }
+
+    //----------  for the correlated Gauss fit (end)   ---------    
+
+
     // reset the display to the starting values
     fDisplay.SetUsed(origPixelsUsed);
@@ -873 +1247 @@
         Float_t starXpos = star->GetMeanX();
         Float_t starYpos = star->GetMeanY();
-        
-        Float_t starSize = 3*star->GetSigmaMajorAxis();
+
+        Float_t starSize = 3*star->GetSigmaMajorAxis();      
         
         Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+
@@ -899 +1273 @@
 
 
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/mtemp/MFindStars.h

@@ -41 +41 @@
     MStarLocalCam *fStars;
 
-    //    MAstroCamera fAstro;
+    MAstroCamera fAstro;
     TArrayC      fPixelsUsed;
     MHCamera     fDisplay;
@@ -54 +54 @@
 
     //Fitting(Minuit) variables
-    const Int_t fNumVar;
+    Int_t fNumVar;
     Float_t fTempChisquare;
     Int_t fTempDegreesofFreedom;
     Int_t fMinuitPrintOutLevel;
     
+    Bool_t fUseCorrelatedGauss;
+    Bool_t fFixCorrelation;
+
     TString *fVname;
     TArrayD fVinit; 
@@ -74 +77 @@
     Bool_t ShadowStar(MStarLocalPos* star);
 
+    TString fGeometryFile;
+    TString fBSCFile;
+
   public:
     
@@ -82 +88 @@
     Int_t PostProcess();
 
+    void SetUseCorrelatedGauss(Bool_t usecorrgauss = kTRUE,
+                               Bool_t fixcorr      = kFALSE);
+    Bool_t  GetUseCorrelatedGauss()          {return fUseCorrelatedGauss;}
+    Bool_t  GetFixCorrelation()              {return fFixCorrelation;}
+
     // setters
-    void SetNumIntegratedEvents(UInt_t max) {fMaxNumIntegratedEvents=max;}
-    void SetRingInterest(Float_t ring) {fRingInterest=ring;}
+    void SetNumIntegratedEvents(UInt_t max)  {fMaxNumIntegratedEvents=max;}
+    void SetRingInterest(Float_t ring)       {fRingInterest=ring;}
     void SetBlindPixels(TArrayS blindpixels);
     void SetMinuitPrintOutLevel(Int_t level) {fMinuitPrintOutLevel=level;}
-    void SetDCTailCut(Float_t cut) {fDCTailCut=cut;}
+    void SetDCTailCut(Float_t cut)           {fDCTailCut=cut;}
 
-    void SetChisquare(Float_t chi) {fTempChisquare=chi;}
-    void SetDegreesofFreedom(Int_t free) {fTempDegreesofFreedom=free;}
+    void SetChisquare(Float_t chi)           {fTempChisquare=chi;}
+    void SetDegreesofFreedom(Int_t free)     {fTempDegreesofFreedom=free;}
+
+    void SetGeometryFile(TString f)          {fGeometryFile=f;}
+    void SetBSCFile(TString f)               {fBSCFile=f;}
+    void SetRaDec(Double_t ra, Double_t dec) {fAstro.SetRaDec(ra,dec);}
+    void SetRaDec(TVector3 &v)               { fAstro.SetRaDec(v); }
+    void SetLimMag(Double_t mag)             { fAstro.SetLimMag(mag); } 
+    void SetRadiusFOV(Double_t deg)          { fAstro.SetRadiusFOV(deg); }
+
 
     //Getters
-    MHCamera& GetDisplay() { return fDisplay; }
+    MHCamera& GetDisplay()           { return fDisplay; }
     
-    Float_t GetChisquare() {return fTempChisquare;}
-    Int_t GetDegreesofFreedom() {return fTempDegreesofFreedom;}
-    UInt_t GetNumIntegratedEvents() {return fMaxNumIntegratedEvents;}
-    Float_t GetRingInterest() {return fRingInterest;}
+    Float_t GetChisquare()           {return fTempChisquare;}
+    Int_t   GetDegreesofFreedom()    {return fTempDegreesofFreedom;}
+    UInt_t  GetNumIntegratedEvents() {return fMaxNumIntegratedEvents;}
+    Float_t GetRingInterest()        {return fRingInterest;}
     
   ClassDef(MFindStars, 0) // Tool to find stars from DC Currents
@@ -104 +123 @@
 
 #endif
+
+