Changeset 4060

Timestamp:

05/13/04 14:03:32 (21 years ago)

Author:

jlopez

Message:

*** empty log message ***

File:

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

trunk/MagicSoft/Mars/mtemp/MFindStars.cc

@@ -43 +43 @@
 #include "MMcConfigRunHeader.h"
 
-#include "MMinuitInterface.h"
-
 #include "MLog.h"
 #include "MLogManip.h"
@@ -84 +82 @@
 
     MFindStars* find =  (MFindStars*)gMinuit->GetObjectFit();
-    MHCamera* display = (MHCamera*)find->GetDisplay();
+    MHCamera& display = (MHCamera&)find->GetDisplay();
     Float_t ped = find->GetPedestalDC();
     Float_t rms = find->GetPedestalRMSDC();
-    MGeomCam& geom = (MGeomCam&)display->GetGeomCam();
+    MGeomCam& geom = (MGeomCam&)display.GetGeomCam();
 
     UInt_t numPixels = geom.GetNumPixels();
@@ -100 +98 @@
     for (UInt_t pixid=1; pixid<numPixels; pixid++) 
     {
-        z = display->GetBinContent(pixid+1)-ped;
-
-        if (display->IsUsed(pixid) && z > 0.)
+        z = display.GetBinContent(pixid+1)-ped;
+
+        if (display.IsUsed(pixid))
         {
             x = geom[pixid].GetX();
@@ -131 +129 @@
   fMaxNumIntegratedEvents = 10;
   fRingInterest = 125.; //[mm] ~ 0.4 deg
-  fMinDCForStars = 2.*fMaxNumIntegratedEvents; //[uA]
+  fDCTailCut = 4;
   
   fPixelsUsed.Set(577);
@@ -138 +136 @@
   //Fitting(Minuit) initialitation
   const Float_t pixelSize = 31.5; //[mm]
-
+  fMinuitPrintOutLevel = -1;
+  
   fVname = new TString[fNumVar];
   fVinit.Set(fNumVar); 
@@ -155 +154 @@
   fVname[1] = "meanx";
   fVinit[1] = 0.;
-  fStep[1]  = fVinit[0]/sqrt2;
+  fStep[1]  = fVinit[1]/sqrt2;
   fLimlo[1] = -600.;
   fLimup[1] = 600.;
@@ -162 +161 @@
   fVname[2] = "sigmaminor";
   fVinit[2] = pixelSize;
-  fStep[2]  = fVinit[0]/sqrt2;
+  fStep[2]  = fVinit[2]/sqrt2;
   fLimlo[2] = pixelSize/(2*sqrt3);
   fLimup[2] = 500.;
@@ -169 +168 @@
   fVname[3] = "meany";
   fVinit[3] = 0.;
-  fStep[3]  = fVinit[0]/sqrt2;
+  fStep[3]  = fVinit[3]/sqrt2;
   fLimlo[3] = -600.;
   fLimup[3] = 600.;
@@ -176 +175 @@
   fVname[4] = "sigmamajor";
   fVinit[4] = pixelSize;
-  fStep[4]  = fVinit[0]/sqrt2;
+  fStep[4]  = fVinit[4]/sqrt2;
   fLimlo[4] = pixelSize/(2*sqrt3);
   fLimup[4] = 500.;
@@ -183 +182 @@
   fObjectFit  = NULL;
   //  fMethod     = "SIMPLEX";
-  //  fMethod     = "MIGRAD";
-  fMethod     = "MINIMIZE";
+  fMethod     = "MIGRAD";
+  //  fMethod     = "MINIMIZE";
   fNulloutput = kFALSE;
 }
@@ -199 +198 @@
     }
 
+    // Initialize camera display with the MGeomCam information
     fDisplay.SetGeometry(*fGeomCam);
     fDisplay.SetUsed(fPixelsUsed);
@@ -281 +281 @@
       return kFALSE;
     }
+
+    fMinDCForStars = 1.5*fMaxNumIntegratedEvents; //[uA]
+
+    // Initialize the TMinuit object
+
+    TMinuit *gMinuit = new TMinuit(fNumVar);  //initialize TMinuit with a maximum of params
+    gMinuit->SetFCN(fcn);
+
+    Double_t arglist[10];
+    Int_t ierflg = 0;
+
+    arglist[0] = 1;
+    gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
+    arglist[0] = fMinuitPrintOutLevel;
+    gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
+
+    // Set object pointer to allow mimuit to access internally
+    gMinuit->SetObjectFit(this);
 
     return kTRUE;
@@ -319 +337 @@
           
           DCPedestalCalc(fPedestalDC, fPedestalRMSDC);
-          fMinDCForStars = fMinDCForStars>(fPedestalDC+5*fPedestalRMSDC)?fMinDCForStars:(fPedestalDC+5*fPedestalRMSDC);
+          fMinDCForStars = fMinDCForStars>(fPedestalDC+fDCTailCut*fPedestalRMSDC)?fMinDCForStars:(fPedestalDC+fDCTailCut*fPedestalRMSDC);
 
           *fLog << dbg << " DC pedestal = " << fPedestalDC << " pedestal rms = " << fPedestalRMSDC << endl;
@@ -416 +434 @@
 Bool_t MFindStars::DCPedestalCalc(Float_t &ped, Float_t &rms)
 {
-  //-------------------------------------------------------------
-  // save pointer to the MINUIT object for optimizing the supercuts
-  // because it will be overwritten
-  // when fitting the alpha distribution in MHFindSignificance
-  TMinuit *savePointer = gMinuit;
-  //-------------------------------------------------------------
+    //-------------------------------------------------------------
+    // save pointer to the MINUIT object for optimizing the supercuts
+    // because it will be overwritten 
+    // when fitting the alpha distribution in MHFindSignificance
+    TMinuit *savePointer = gMinuit;
+    //-------------------------------------------------------------
 
    UInt_t numPixels = fGeomCam->GetNumPixels();
@@ -453 +471 @@
    
    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());
@@ -465 +486 @@
    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;
-   //-------------------------------------------
+    //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;
+    //-------------------------------------------
 
    return kTRUE;
@@ -598 +619 @@
             usedPx++;
 
-            Float_t charge  = fDisplay.GetBinContent(pix+1);
+            Float_t charge  = fDisplay.GetBinContent(pix+1)-fPedestalDC;
             Float_t pixXpos = (*fGeomCam)[pix].GetX();
             Float_t pixYpos = (*fGeomCam)[pix].GetY();
@@ -618 +639 @@
     meanSqY /= sumCharge;
 
-    //Substract pedestal from sumCharge
-    sumCharge -= usedPx*fPedestalDC;
-
-    Float_t rmsX = TMath::Sqrt(meanSqX - meanX*meanX);
-    Float_t rmsY = TMath::Sqrt(meanSqY - meanY*meanY);
+    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 << " 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 << endl;
+        rmsY = 0.;
+      }
     
     star->SetCalcValues(sumCharge,max,meanX,meanY,rmsX,rmsY);
 
-
+    if (rmsX <= 0. || rmsY <= 0.)
+      return kFALSE;
+    
 // fit the star spot using TMinuit
 
-  for (UInt_t pix=1; pix<numPixels; pix++)
+    
+    for (UInt_t pix=1; pix<numPixels; pix++)
       if (fDisplay.IsUsed(pix))
         *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl;
   
   //Initialate variables for fit
-    fVinit[0] = star->GetMaxCalc()-fPedestalDC;
+    fVinit[0] = max;
     fLimlo[0] = fMinDCForStars-fPedestalDC;
-    fLimup[0] = fLimup[0]-fPedestalDC;
+    fLimup[0] = 30*fMaxNumIntegratedEvents-fPedestalDC;
     fVinit[1] = meanX;
     fVinit[2] = rmsX;
@@ -647 +686 @@
     //
 
+    // -------------------------------------------
+    // call MINUIT
+
+    Double_t arglist[10];
+    Int_t ierflg = 0;
+
+    for (Int_t i=0; i<fNumVar; i++)
+      gMinuit->mnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg);
+
     TStopwatch clock;
     clock.Start();
 
-    *fLog << dbg << " before calling CallMinuit" << endl;
-
-    MMinuitInterface inter;               
-    Bool_t rc = inter.CallMinuit(fcn, fVname,
-                                 fVinit, fStep, fLimlo, fLimup, fFix,
-                                 this, fMethod, fNulloutput);
- 
-    *fLog << dbg << "after calling CallMinuit" << endl;
-    *fLog << dbg << "Time spent for the minimization in MINUIT :   " << endl;;
+// Now ready for minimization step
+    arglist[0] = 500;
+    arglist[1] = 1.;
+    gMinuit->mnexcm(fMethod, arglist ,2,ierflg);
+
     clock.Stop();
-    clock.Print();
+
+    if(fMinuitPrintOutLevel>=0)
+      {
+        *fLog << dbg << "Time spent for the minimization in MINUIT :   " << endl;;
+        clock.Print();
+      }
 
     Double_t integratedCharge;
@@ -671 +720 @@
     Int_t   dregrees  = GetDegreesofFreedom()-fNumVar;
 
-    if (rc)
+    if (!ierflg)
       {
         gMinuit->GetParameter(0,maxFit, maxFitError);
@@ -692 +741 @@
         sigmaMajorAxis = 0.;
         integratedCharge = 0.;
-      }
-    
-    
+
+        *fLog << err << "TMinuit::Call error " << ierflg << endl;
+      }
     
     star->SetFitValues(integratedCharge,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,chisquare,dregrees);
-
+    
     // reset the display to the starting values
     fDisplay.SetUsed(origPixelsUsed);
 
-    return rc;
+    if (ierflg)
+      return kCONTINUE;
+    return kTRUE;
 }
 
@@ -740 +791 @@
 
 
+
+