Changeset 4294

Timestamp:

06/16/04 09:57:59 (20 years ago)

Author:

jlopez

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp

Files:

• MFindStars.cc (modified) (33 diffs)
• MFindStars.h (modified) (5 diffs)
• MStarLocalCam.cc (modified) (3 diffs)
• MStarLocalCam.h (modified) (3 diffs)
• MStarLocalPos.cc (modified) (1 diff)
• MStarLocalPos.h (modified) (1 diff)
• mifae/Changelog (modified) (1 diff)
• mifae/library/IFAELinkDef.h (modified) (1 diff)
• mifae/library/MCalibrateDC.cc (modified) (1 diff)
• mifae/library/Makefile (modified) (2 diffs)
• mifae/macros/findstars.C (modified) (2 diffs)

trunk/MagicSoft/Mars/mtemp/MFindStars.cc

@@ -36 +36 @@
 #include <TFile.h>
 #include <TTree.h>
+#include <TCanvas.h>
 #include <TH1F.h>
 #include <TF1.h>
@@ -57 +58 @@
 
 #include "MParList.h"
+#include "MTaskList.h"
 
 ClassImp(MFindStars);
@@ -63 +65 @@
 const Float_t sqrt2 = sqrt(2.);
 const Float_t sqrt3 = sqrt(3.);
+const UInt_t  lastInnerPixel = 396;
+    
 
 //______________________________________________________________________________
@@ -81 +85 @@
 {
 
-    MFindStars* find =  (MFindStars*)gMinuit->GetObjectFit();
-    MHCamera& display = (MHCamera&)find->GetDisplay();
-    Float_t ped = find->GetPedestalDC();
-    Float_t rms = find->GetPedestalRMSDC();
-    MGeomCam& geom = (MGeomCam&)display.GetGeomCam();
-
-    UInt_t numPixels = geom.GetNumPixels();
- 
+  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");
+
+  MHCamera& display = (MHCamera&)find->GetDisplay();
+  
+  Float_t innerped = stars->GetInnerPedestalDC();
+  Float_t innerrms = stars->GetInnerPedestalRMSDC();
+  Float_t outerped = stars->GetOuterPedestalDC();
+  Float_t outerrms = stars->GetOuterPedestalRMSDC();
+
+  UInt_t numPixels = geom->GetNumPixels();
+  
 //calculate chisquare
     Double_t chisq = 0;
     Double_t delta;
     Double_t x,y,z;
-    Double_t errorz = rms; //[uA]
+    Double_t errorz=0;
 
     UInt_t usedPx=0;
     for (UInt_t pixid=1; pixid<numPixels; pixid++) 
     {
-        z = display.GetBinContent(pixid+1)-ped;
-
         if (display.IsUsed(pixid))
         {
-            x = geom[pixid].GetX();
-            y = geom[pixid].GetY();
+            x = (*geom)[pixid].GetX();
+            y = (*geom)[pixid].GetY();
+            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;
+              delta  = (z-func(x,y,par))/errorz;
+              chisq += delta*delta;
             }
             else
@@ -121 +132 @@
 }
 
-MFindStars::MFindStars(const char *name, const char *title): fNumVar(5)
+MFindStars::MFindStars(const char *name, const char *title): 
+  fGeomCam(NULL), fCurr(NULL), fTimeCurr(NULL), fDrive(NULL), fStars(NULL), fNumVar(5)
 {
   fName  = name  ? name  : "MFindStars";
@@ -185 +197 @@
   //  fMethod     = "MINIMIZE";
   fNulloutput = kFALSE;
+
+  // Set output level
+  //  fLog->SetOutputLevel(3); // No dbg messages
 }
 
@@ -199 +214 @@
 
     // Initialize camera display with the MGeomCam information
-    fDisplay.SetGeometry(*fGeomCam);
+    fDisplay.SetGeometry(*fGeomCam,"FindStarsDisplay","FindStarsDisplay");
     fDisplay.SetUsed(fPixelsUsed);
-    
+
     fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));
 
@@ -218 +233 @@
     }
 
-    fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
+    //    fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive"));
 
     if (!fDrive)
@@ -226 +241 @@
     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());
+//         //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());
         
       }
@@ -282 +297 @@
     }
 
-    fMinDCForStars = 1.5*fMaxNumIntegratedEvents; //[uA]
+    fMinDCForStars = 1.*fMaxNumIntegratedEvents; //[uA]
 
     // Initialize the TMinuit object
@@ -297 +312 @@
     gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
 
-    // Set object pointer to allow mimuit to access internally
-    gMinuit->SetObjectFit(this);
+    // Set MParList object pointer to allow mimuit to access internally
+    gMinuit->SetObjectFit(pList);
 
     return kTRUE;
@@ -305 +320 @@
 Int_t MFindStars::Process()
 {
+
   UInt_t numPixels = fGeomCam->GetNumPixels();
   TArrayC origPixelsUsed;
@@ -314 +330 @@
       if (fDrive)
         {
-          Float_t ra  = fDrive->GetRa();
-          Float_t dec = fDrive->GetDec();
+//           Float_t ra  = fDrive->GetRa();
+//           Float_t dec = fDrive->GetDec();
           
-          fAstro.SetRaDec(ra, dec);
-          fAstro.SetGuiActive();
+//           fAstro.SetRaDec(ra, dec);
+//           fAstro.SetGuiActive();
           
-          fAstro.FillStarList();
+//           fAstro.FillStarList();
         }
       else
@@ -336 +352 @@
             }
           
-          DCPedestalCalc(fPedestalDC, fPedestalRMSDC);
-          fMinDCForStars = fMinDCForStars>(fPedestalDC+fDCTailCut*fPedestalRMSDC)?fMinDCForStars:(fPedestalDC+fDCTailCut*fPedestalRMSDC);
-
-          *fLog << dbg << " DC pedestal = " << fPedestalDC << " pedestal rms = " << fPedestalRMSDC << endl;
-          *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))
+          if (DCPedestalCalc())
             {
 
-              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);
+              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);
             }
-
-          fDisplay.SetUsed(origPixelsUsed);
+          
         }
       
@@ -366 +386 @@
         {
           *fLog << err << GetName() << " No stars candidates in the camera." << endl;
-          return kFALSE;
+          return kCONTINUE;
         }
       else
@@ -390 +410 @@
       //After finding stars reset all vairables
       fDisplay.Reset();
+      fStars->GetDisplay().Reset(); //FIXME: Put this display just in the container
       fDisplay.SetUsed(origPixelsUsed);
       fNumIntegratedEvents=0;
@@ -399 +420 @@
           origPixelsUsed[pix]=(Char_t)kTRUE;
         else
-            origPixelsUsed[pix]=(Char_t)kFALSE;
+          origPixelsUsed[pix]=(Char_t)kFALSE;
         
       }
@@ -413 +434 @@
 Int_t MFindStars::PostProcess()
 {
-  if(fStars->GetList()->GetSize() != 0)
-    fStars->Print();
-
   return kTRUE;
 }
@@ -426 +444 @@
       {
         fPixelsUsed[blindpixels[idx]]=(Char_t)kFALSE;
-        *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <<blindpixels[idx]  << ") kFALSE" << endl;
+        if (fMinuitPrintOutLevel>=0) *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <<blindpixels[idx]  << ") kFALSE" << endl;
       }
     
@@ -432 +450 @@
 }
 
-Bool_t MFindStars::DCPedestalCalc(Float_t &ped, Float_t &rms)
+Bool_t MFindStars::DCPedestalCalc()
 {
     //-------------------------------------------------------------
@@ -442 +460 @@
 
    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);
+   Float_t ped;
+   Float_t rms;
+
+   TH1F **dchist = new TH1F*[2];
+   for (UInt_t i=0; i<2; i++)
+      dchist[i] = new TH1F("","",26,0.4*fMaxNumIntegratedEvents,3.*fMaxNumIntegratedEvents);
    
-   dchist.Fit("func","QR0");
-   // Remove the comments if you want to go through the file
-   // event-by-event:
-   //   HandleInput();
-
-   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;
-
-    //Check results from the fit are consistent
-    if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess)
-      {
-        *fLog << warn << GetName() << " Pedestal DC fit give non consistent results." << endl;
-        *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
-        *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
-        ped = maxprobdc;
-        rms = rmsguess/1.175; // FWHM=2.35*rms
-      }
-
-    //=================================================================
-
-    // reset gMinuit to the MINUIT object for optimizing the supercuts 
-    gMinuit = savePointer;
-    //-------------------------------------------
-
+   for (UInt_t pix=1; pix<=lastInnerPixel; pix++)
+       dchist[0]->Fill(fDisplay.GetBinContent(pix+1));
+   for (UInt_t pix=lastInnerPixel+1; pix<numPixels; pix++)
+       dchist[1]->Fill(fDisplay.GetBinContent(pix+1));
+
+   // 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());
+      UInt_t bin = dchist[i]->GetMaximumBin();
+      do
+        {
+          bin++;
+        }
+      while(dchist[i]->GetBinContent(bin)/nummaxprobdc > 0.5);
+      Float_t halfmaxprobdc = dchist[i]->GetBinCenter(bin);
+      
+      if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
+      if (fMinuitPrintOutLevel>=0) *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist[i]->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;
+      
+      if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
+      
+      TF1 func("func","gaus",min,max);
+      func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
+      
+      dchist[i]->Fit("func","QR0");
+      
+      UInt_t aproxnumdegrees = 6*(bin-dchist[i]->GetMaximumBin());
+      Float_t chiq = func.GetChisquare();
+      ped = func.GetParameter(1);
+      rms = func.GetParameter(2);
+      
+      if (fMinuitPrintOutLevel>=0) *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
+      
+      Int_t numsigmas = 5;
+      Axis_t minbin = ped-numsigmas*rms/dchist[i]->GetBinWidth(1);
+      minbin=minbin<1?1:minbin;
+      Axis_t maxbin = ped+numsigmas*rms/dchist[i]->GetBinWidth(1);
+      if (fMinuitPrintOutLevel>=0) *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist[i]->Integral((int)minbin,(int)maxbin) << endl;
+      
+      //Check results from the fit are consistent
+      if (ped < 0. || rms < 0.)
+        {
+          *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)
+        {
+          *fLog << warn << GetName() << " Pedestal DC fit give non consistent results for " << (i==0?"Inner":"Outer") << "pixels." << endl;
+          *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
+          *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
+          ped = maxprobdc;
+          rms = rmsguess/1.175; // FWHM=2.35*rms
+        }
+   
+      if (i == 0)
+        {
+          fStars->SetInnerPedestalDC(ped);
+          fStars->SetInnerPedestalRMSDC(rms);
+        }
+      else
+        {
+          fStars->SetOuterPedestalDC(ped);
+          fStars->SetOuterPedestalRMSDC(rms);
+        }
+
+      
+
+    }
+   
+
+   for (UInt_t i=0; i<2; i++)
+      delete dchist[i];
+   delete [] dchist;
+
+   //=================================================================
+
+   // reset gMinuit to the MINUIT object for optimizing the supercuts 
+   gMinuit = savePointer;
+   //-------------------------------------------
+   
+   if (fStars->GetInnerPedestalDC() < 0. ||  fStars->GetInnerPedestalRMSDC() < 0. || fStars->GetOuterPedestalDC() < 0. ||  fStars->GetOuterPedestalRMSDC() < 0.)
+     return kFALSE;
+  
    return kTRUE;
 }
@@ -561 +618 @@
 
     if (maxDC == 0)
+      {
+        *fLog << warn << " No found pixels with maximum dc" << endl;
         return kFALSE;
-
+      }
+    
     maxPix = (*fGeomCam)[maxPixIdx[0]];
 
-    *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxDC << " uA) x position(" << setw(3) << maxPix.GetX() <<  " mm) x position(" << setw(3) << maxPix.GetY() << " mm) swnumber(" << maxPixIdx[0] << ")" << endl;
+    if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxDC << " uA) x position(" << setw(3) << maxPix.GetX() <<  " mm) x position(" << setw(3) << maxPix.GetY() << " mm) swnumber(" << maxPixIdx[0] << ")" << endl;
 
     return kTRUE;
@@ -574 +634 @@
 
   UInt_t numPixels = fGeomCam->GetNumPixels();
+  Float_t innerped = fStars->GetInnerPedestalDC();
+  Float_t outerped = fStars->GetOuterPedestalDC();
+
   TArrayC origPixelsUsed;
   origPixelsUsed.Set(numPixels);
@@ -589 +652 @@
   
   Float_t max=0;
+  UInt_t  pixmax=0;
   Float_t meanX=0;
   Float_t meanY=0;
@@ -594 +658 @@
   Float_t meanSqY=0;
   Float_t sumCharge=0;
-  UInt_t usedPx=0;      
+  UInt_t usedInnerPx=0; 
+  UInt_t usedOuterPx=0; 
 
   const Float_t meanPresition = 3.; //[mm]
@@ -620 +685 @@
     
 // determine mean x and mean y
-    usedPx=0;   
+    usedInnerPx=0;      
+    usedOuterPx=0;      
     for(UInt_t pix=0; pix<numPixels; pix++)
     {
         if(fDisplay.IsUsed(pix))
         {
-            usedPx++;
+            pix>lastInnerPixel?usedOuterPx++:usedInnerPx++;
 
             Float_t charge  = fDisplay.GetBinContent(pix+1);
@@ -631 +697 @@
             Float_t pixYpos = (*fGeomCam)[pix].GetY();
 
-            if (charge>max) max=charge;
+            if (charge>max)
+              {
+                max=charge;
+                pixmax=pix;
+              }
+            
             meanX     += charge*pixXpos;
             meanY     += charge*pixYpos;
@@ -640 +711 @@
     }
 
-    *fLog << dbg << " usedPx " << usedPx << endl;
+    if (fMinuitPrintOutLevel>=0) *fLog << dbg << " usedInnerPx " << usedInnerPx << " usedOuterPx " << usedOuterPx << endl;
 
     meanX   /= sumCharge;
@@ -680 +751 @@
     
     // Substrack pedestal DC
-    sumCharge-=fPedestalDC*usedPx;
-    max-=fPedestalDC;
+    sumCharge-= (usedInnerPx*innerped+usedOuterPx*outerped)/(usedInnerPx+usedOuterPx);
+    max-=pixmax>lastInnerPixel?outerped:innerped;
     
 
@@ -695 +766 @@
     for (UInt_t pix=1; pix<numPixels; pix++)
       if (fDisplay.IsUsed(pix))
-        *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
-  
+        if (fMinuitPrintOutLevel>=0) *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
+  
+    if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars " << fMinDCForStars << " pixmax>lastInnerPixel?outerped:innerped " << (pixmax>lastInnerPixel?outerped:innerped) << " fMaxNumIntegratedEvents " << fMaxNumIntegratedEvents << endl;
+
   //Initialate variables for fit
     fVinit[0] = max;
-    fLimlo[0] = fMinDCForStars-fPedestalDC;
-    fLimup[0] = 30*fMaxNumIntegratedEvents-fPedestalDC;
+    fLimlo[0] = fMinDCForStars-(pixmax>lastInnerPixel?outerped:innerped);
+    fLimup[0] = 30*fMaxNumIntegratedEvents-(pixmax>lastInnerPixel?outerped:innerped);
     fVinit[1] = meanX;
     fVinit[2] = rmsX;
@@ -707 +780 @@
     //Init steps
     for(Int_t i=0; i<fNumVar; i++)
+      {
         if (fVinit[i] != 0)
           fStep[i] = TMath::Abs(fVinit[i]/sqrt2);
+        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;
+      }
     //
 
@@ -732 +811 @@
     if(fMinuitPrintOutLevel>=0)
       {
-        *fLog << dbg << "Time spent for the minimization in MINUIT :   " << endl;;
+        if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Time spent for the minimization in MINUIT :   " << endl;;
         clock.Print();
       }
@@ -807 +886 @@
     }
 
-    *fLog << dbg << " shadowPx " << shadowPx << endl;
+    if (fMinuitPrintOutLevel>=0) *fLog << dbg << " shadowPx " << shadowPx << endl;
 
     fDisplay.SetUsed(fPixelsUsed);

trunk/MagicSoft/Mars/mtemp/MFindStars.h

@@ -41 +41 @@
     MStarLocalCam *fStars;
 
-    MAstroCamera fAstro;
+    //    MAstroCamera fAstro;
     TArrayC      fPixelsUsed;
     MHCamera     fDisplay;
@@ -51 +51 @@
     Float_t fMinDCForStars; //[uA]
 
-    Float_t fPedestalDC;    //[ua]
-    Float_t fPedestalRMSDC; //[ua]
+    Float_t fDCTailCut;
 
     //Fitting(Minuit) variables
@@ -58 +57 @@
     Float_t fTempChisquare;
     Int_t fTempDegreesofFreedom;
+    Int_t fMinuitPrintOutLevel;
     
     TString *fVname;
@@ -69 +69 @@
     Bool_t fNulloutput;
     
-    Bool_t DCPedestalCalc(Float_t &ped, Float_t &rms);
+    Bool_t DCPedestalCalc();
     Bool_t FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix);
     Bool_t FindStar(MStarLocalPos* star);
@@ -81 +81 @@
     Int_t Process();
     Int_t PostProcess();
-  
+
+    // setters
     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 SetChisquare(Float_t chi) {fTempChisquare=chi;}
     void SetDegreesofFreedom(Int_t free) {fTempDegreesofFreedom=free;}
 
-    MHCamera* GetDisplay() { return &fDisplay; }
-    Float_t GetPedestalDC() { return fPedestalDC; }
-    Float_t GetPedestalRMSDC() { return fPedestalRMSDC; }
+    //Getters
+    MHCamera& GetDisplay() { return fDisplay; }
     
     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

trunk/MagicSoft/Mars/mtemp/MStarLocalCam.cc

@@ -53 +53 @@
   
   fStars = new TList;
+
+  fInnerPedestalDC = 0.;
+  fOuterPedestalDC = 0.;
+  fInnerPedestalRMSDC = 0.;
+  fOuterPedestalRMSDC = 0.;
+  
 }
 
@@ -67 +73 @@
 MStarLocalPos &MStarLocalCam::operator[] (Int_t i)
 {
-    return *static_cast<MStarLocalPos*>(fStars->At(i));
+  MStarLocalPos& star = *static_cast<MStarLocalPos*>(fStars->At(i));
+  return star;
 }
 
@@ -89 +96 @@
 void MStarLocalCam::Print(Option_t *o) const
 {
-        TIter Next(fStars);
-        MStarLocalPos* star;
-        UInt_t starnum = 0;
-        while ((star=(MStarLocalPos*)Next())) 
-          {
-            *fLog << inf << "Star[" << starnum << "] info:" << endl;
-            star->Print(o);
-            starnum++;
-          }
+      *fLog << inf << "DCs baseline:" << endl;
+      *fLog << inf << " Inner Pixels Mean pedestal \t" << setw(4) << fInnerPedestalDC << " uA  and RMS " << setw(4) << fInnerPedestalRMSDC << " uA for DCs" << endl;
+      *fLog << inf << " Outer Pixels Mean pedestal \t" << setw(4) << fOuterPedestalDC << " uA  and RMS " << setw(4) << fOuterPedestalRMSDC << " uA for DCs" << endl;
+      
+      TIter Next(fStars);
+      MStarLocalPos* star;
+      UInt_t starnum = 0;
+      while ((star=(MStarLocalPos*)Next())) 
+        {
+          *fLog << inf << "Star[" << starnum << "] info:" << endl;
+          star->Print(o);
+          starnum++;
+        }
 }

trunk/MagicSoft/Mars/mtemp/MStarLocalCam.h

@@ -5 +5 @@
 #include "MParContainer.h"
 #endif
-#ifndef MARS_MCamEvent
-#include "MCamEvent.h"
+
+#ifndef MARS_MHCamera
+#include "MHCamera.h"
 #endif
 
@@ -20 +21 @@
   TList  *fStars;  //-> FIXME: Change TClonesArray away from a pointer?
 
-public:
+  // BaseLine DCs info
+  Float_t fInnerPedestalDC;         //[ua]
+  Float_t fOuterPedestalDC;         //[ua]
+
+  Float_t fInnerPedestalRMSDC;      //[ua]
+  Float_t fOuterPedestalRMSDC;      //[ua]
+
+  MHCamera fDisplay;
+
+ public:
 
   MStarLocalCam(const char *name=NULL, const char *title=NULL);
@@ -29 +39 @@
 
   TList *GetList() const { return fStars; }
+  UInt_t GetNumStars() const { return fStars->GetSize(); }
+
+  //Getters
+
+  Float_t GetInnerPedestalDC() {return fInnerPedestalDC;}
+  Float_t GetOuterPedestalDC() {return fOuterPedestalDC;}
+  Float_t GetInnerPedestalRMSDC() { return fInnerPedestalRMSDC;}
+  Float_t GetOuterPedestalRMSDC() { return fOuterPedestalRMSDC;}
+
+  MHCamera& GetDisplay() { return fDisplay; }
+
+    //Setters
+  void SetInnerPedestalDC(Float_t ped) {fInnerPedestalDC = ped;}
+  void SetOuterPedestalDC(Float_t ped) {fOuterPedestalDC = ped;}
+  void SetInnerPedestalRMSDC(Float_t rms){fInnerPedestalRMSDC = rms;}
+  void SetOuterPedestalRMSDC(Float_t rms){fOuterPedestalRMSDC = rms;}
 
   void Paint(Option_t *o=NULL);

trunk/MagicSoft/Mars/mtemp/MStarLocalPos.cc

@@ -68 +68 @@
      fSigmaMajorAxisFit = 0.;
      fChiSquare = 0.;
-     fNdof = 0;
+     fNdof = 1;
 
 }

trunk/MagicSoft/Mars/mtemp/MStarLocalPos.h

@@ -59 +59 @@
     Float_t GetSigmaMajorAxisFit() {return fSigmaMajorAxisFit;}
     Float_t GetChiSquare() {return fChiSquare;}
-    Float_t GetNdof() {return fNdof;}
+    UInt_t GetNdof() {return fNdof;}
     Float_t GetChiSquareNdof() {return fChiSquare/fNdof;}
 

trunk/MagicSoft/Mars/mtemp/mifae/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+  2004/06/15  Javier Lopez
+    * library/MSrcPosFromStars.[h,cc]
+      - Task to compute the position of a source from
+        the positions of several stars.
+
+    * library/MHPSFFromStars.[h,cc]
+      - Histogram task that holds and fills the histograms of 
+        positioning and point spread function of an star.
+
+    * macros/distancebetweenstars.C
+      - Add macros that show in an histogram the distance between
+        the stars in a field of view. 
+
+    * macros/psffromstars.C
+      - Add task that show the point spread funtion and positon of
+        the most central star in the field of view. This information
+        is compute using the MFindStars tool.
 
  2004/06/11 Ester Aliu 

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

@@ -18 +18 @@
 #pragma link C++ class MFHVNotNominal+;
 #pragma link C++ class MCalibrateDC+;
+#pragma link C++ class MHPSFFromStars+;
+#pragma link C++ class MSrcPosFromStars+;
 
 #endif

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

@@ -51 +51 @@
 using namespace std;
 
-MCalibrateDC::MCalibrateDC(TString filename, const char *name, const char *title) : fMinDCAllowed(0.5) //[ua]
+MCalibrateDC::MCalibrateDC(TString filename, const char *name, const char *title) 
 {
   fName  = name  ? name  : "MCalibrateDC";

trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile

@@ -44 +44 @@
            -I../../../mcamera \
            -I../../../mastro \
-           -I../../../mhist
+           -I../../../mhist \
+           -I../../
 
 
@@ -60 +61 @@
         MIslandClean.cc \
         MFHVNotNominal.cc \
-        MCalibrateDC.cc
+        MCalibrateDC.cc \
+        MHPSFFromStars.cc \
+        MSrcPosFromStars.cc
 
 ############################################################

trunk/MagicSoft/Mars/mtemp/mifae/macros/findstars.C

@@ -47 +47 @@
 
 
-void findstars(const TString filename="dc_2004_03_19_00_36_50_20781_Mrk421.root", const TString directory="/nfs/magic/CaCodata/rootdata/Mrk421/Period015/2004_03_19/", const UInt_t numEvents = 0)
+void findstars(const TString filename="dc_2004_03_17_01_16_51_20440_Mrk421.root", const TString directory="/nfs/magic/CaCodata/rootdata/Mrk421/Period015/2004_03_17/", const UInt_t numEvents = 0)
 {
 
@@ -81 +81 @@
   MGeomApply geomapl;
   TString continuoslightfile = 
-    "/nfs/magic/CaCodata/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
+    //    "/home/Javi/mnt_magic_data/CaCo/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
+     "/nfs/magic/CaCodata/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
+
+  Float_t mindc = 0.9; //[uA]
   MCalibrateDC dccal;
   dccal.SetFileName(continuoslightfile);
-  
+  dccal.SetMinDCAllowed(mindc);
 
-  const Int_t numblind = 12;
-  const Short_t x[numblind] = {  8,  27,
-                               507, 508, 509, 510, 511,
-                               543, 559, 560, 561, 567};
+  const Int_t numblind = 5;
+  const Short_t x[numblind] = { 47, 124, 470, 475, 571};
   const TArrayS blindpixels(numblind,(Short_t*)x);
   Float_t ringinterest = 100; //[mm]
-  Float_t tailcut = 3.5;
-  UInt_t integratedevents = 10;
+  Float_t tailcut = 2.5;
+  UInt_t integratedevents = 1;
 
   MFindStars findstars;
-  //  findstars.SetBlindPixels(blindpixels);
+  findstars.SetBlindPixels(blindpixels);
   findstars.SetRingInterest(ringinterest);
   findstars.SetDCTailCut(tailcut);
   findstars.SetNumIntegratedEvents(integratedevents);
-  findstars.SetMinuitPrintOutLevel(0);
+  findstars.SetMinuitPrintOutLevel(-1);
 
-  // prints
-  MPrint pdc("MCameraDC");
-  MPrint pstar("MStarLocalCam");
   
   tlist.AddToList(&geomapl);
   tlist.AddToList(&read);
   tlist.AddToList(&dccal);
-  //  tlist.AddToList(&pdc, "Currents");
   tlist.AddToList(&findstars, "Currents");
-  //  tlist.AddToList(&pstar, "Currents");
 
   //