Index: trunk/MagicSoft/Mars/mtemp/MFindStars.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MFindStars.cc (revision 4059) +++ trunk/MagicSoft/Mars/mtemp/MFindStars.cc (revision 4060) @@ -43,6 +43,4 @@ #include "MMcConfigRunHeader.h" -#include "MMinuitInterface.h" - #include "MLog.h" #include "MLogManip.h" @@ -84,8 +82,8 @@ 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,7 +98,7 @@ for (UInt_t pixid=1; pixidGetBinContent(pixid+1)-ped; - - if (display->IsUsed(pixid) && z > 0.) + z = display.GetBinContent(pixid+1)-ped; + + if (display.IsUsed(pixid)) { x = geom[pixid].GetX(); @@ -131,5 +129,5 @@ fMaxNumIntegratedEvents = 10; fRingInterest = 125.; //[mm] ~ 0.4 deg - fMinDCForStars = 2.*fMaxNumIntegratedEvents; //[uA] + fDCTailCut = 4; fPixelsUsed.Set(577); @@ -138,5 +136,6 @@ //Fitting(Minuit) initialitation const Float_t pixelSize = 31.5; //[mm] - + fMinuitPrintOutLevel = -1; + fVname = new TString[fNumVar]; fVinit.Set(fNumVar); @@ -155,5 +154,5 @@ fVname[1] = "meanx"; fVinit[1] = 0.; - fStep[1] = fVinit[0]/sqrt2; + fStep[1] = fVinit[1]/sqrt2; fLimlo[1] = -600.; fLimup[1] = 600.; @@ -162,5 +161,5 @@ 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,5 +168,5 @@ fVname[3] = "meany"; fVinit[3] = 0.; - fStep[3] = fVinit[0]/sqrt2; + fStep[3] = fVinit[3]/sqrt2; fLimlo[3] = -600.; fLimup[3] = 600.; @@ -176,5 +175,5 @@ 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,6 +182,6 @@ fObjectFit = NULL; // fMethod = "SIMPLEX"; - // fMethod = "MIGRAD"; - fMethod = "MINIMIZE"; + fMethod = "MIGRAD"; + // fMethod = "MINIMIZE"; fNulloutput = kFALSE; } @@ -199,4 +198,5 @@ } + // Initialize camera display with the MGeomCam information fDisplay.SetGeometry(*fGeomCam); fDisplay.SetUsed(fPixelsUsed); @@ -281,4 +281,22 @@ 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,5 +337,5 @@ 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,10 +434,10 @@ 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,4 +471,7 @@ 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,21 +486,21 @@ 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,5 +619,5 @@ 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,23 +639,41 @@ 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; pixGetMaxCalc()-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,18 +686,28 @@ // + // ------------------------------------------- + // call MINUIT + + Double_t arglist[10]; + Int_t ierflg = 0; + + for (Int_t i=0; imnparm(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,5 +720,5 @@ Int_t dregrees = GetDegreesofFreedom()-fNumVar; - if (rc) + if (!ierflg) { gMinuit->GetParameter(0,maxFit, maxFitError); @@ -692,14 +741,16 @@ 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,2 +791,4 @@ + +