Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 4545) +++ /trunk/MagicSoft/Mars/Changelog (revision 4546) @@ -54,4 +54,7 @@ and MHTelAxisFromStars was added + + * mtemp/MFindStars.[h,cc] + - add correlated Gauss fit Index: /trunk/MagicSoft/Mars/mtemp/MFindStars.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MFindStars.cc (revision 4545) +++ /trunk/MagicSoft/Mars/mtemp/MFindStars.cc (revision 4546) @@ -15,6 +15,7 @@ ! * ! -! Author(s): Robert Wagner, 2/2004 -! Author(s): Javier López , 4/2004 +! Author(s): Robert Wagner, 2/2004 +! Javier López , 4/2004 +! Wolfgang Wittek, 8/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 @@ -39,4 +40,6 @@ #include #include +#include + #include "MObservatory.h" @@ -70,9 +73,10 @@ //______________________________________________________________________________ // -// 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,4 +84,19 @@ //______________________________________________________________________________ // +// 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,7 +106,8 @@ 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,4 +120,17 @@ 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,16 +142,38 @@ for (UInt_t pixid=1; pixidlastInnerPixel?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,11 +185,25 @@ 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,4 +218,5 @@ const Float_t pixelSize = 31.5; //[mm] fMinuitPrintOutLevel = -1; + //fMinuitPrintOutLevel = 3; fVname = new TString[fNumVar]; @@ -171,5 +241,5 @@ fFix[1] = 0; - fVname[2] = "sigmaminor"; + fVname[2] = "sigmax"; fVinit[2] = pixelSize; fStep[2] = fVinit[2]/sqrt2; @@ -185,5 +255,5 @@ fFix[3] = 0; - fVname[4] = "sigmamajor"; + fVname[4] = "sigmay"; fVinit[4] = pixelSize; fStep[4] = fVinit[4]/sqrt2; @@ -192,4 +262,14 @@ 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,5 +280,34 @@ // 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,62 +342,36 @@ } - // 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,4 +395,10 @@ 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,4 +406,12 @@ return kTRUE; } + +//------------------------------------------------------------------------ +// +// Process : +// +// +// +// Int_t MFindStars::Process() @@ -325,108 +422,127 @@ 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 " <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; pixGetInnerPedestalDC()+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; pixGetList()->Delete(); - // - - for (UInt_t pix=1; pixGetInnerPedestalDC()+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; pixGetList()); - 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; pixGetList()); + 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; pixFill(fDisplay.GetBinContent(pix+1)); @@ -474,5 +608,5 @@ // 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,4 +651,9 @@ { *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,11 +666,4 @@ } -// TCanvas *c1 = new TCanvas("c2","c2",500,800); -// dchist[i]->Draw(); -// c1->Print("dchist.ps"); -// getchar(); -// delete c1; -// exit(1); - if (i == 0) { @@ -544,13 +676,13 @@ fStars->SetOuterPedestalRMSDC(rms); } - - - - } + } + for (UInt_t i=0; i<2; i++) + { delete dchist[i]; - delete [] dchist; + } + //delete [] dchist; //================================================================= @@ -632,7 +764,14 @@ } + +//---------------------------------------------------------------------------- +// +// FindStar +// +// +// + Bool_t MFindStars::FindStar(MStarLocalPos* star) { - UInt_t numPixels = fGeomCam->GetNumPixels(); Float_t innerped = fStars->GetInnerPedestalDC(); @@ -652,5 +791,5 @@ Float_t expX = star->GetXExp(); Float_t expY = star->GetYExp(); - + Float_t max=0; UInt_t pixmax=0; @@ -660,6 +799,9 @@ 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,98 +809,166 @@ UInt_t numIterations = 0; - do - { - // First define a area of interest around the expected position of the star - for (UInt_t pix=1; pixmax) + { + max=charge; + pixmax=pix; + } + } + + // 3.) make it new center + + expX = (*fGeomCam)[pixmax].GetX(); + expY = (*fGeomCam)[pixmax].GetY(); + + for (UInt_t pix=1; pixlastInnerPixel?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,5 +982,5 @@ 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,13 +992,23 @@ //Init steps for(Int_t i=0; i=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,4 +1029,9 @@ 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,40 +1042,189 @@ } - 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; kGetParameter( 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; jmnpout(j, name, val, err, bnd1, bnd2, jvarbl); + + for (Int_t k=0; kmnpout(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=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,6 +1247,6 @@ 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,2 +1273,9 @@ + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/MFindStars.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MFindStars.h (revision 4545) +++ /trunk/MagicSoft/Mars/mtemp/MFindStars.h (revision 4546) @@ -41,5 +41,5 @@ MStarLocalCam *fStars; - // MAstroCamera fAstro; + MAstroCamera fAstro; TArrayC fPixelsUsed; MHCamera fDisplay; @@ -54,9 +54,12 @@ //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,4 +77,7 @@ Bool_t ShadowStar(MStarLocalPos* star); + TString fGeometryFile; + TString fBSCFile; + public: @@ -82,21 +88,34 @@ 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,2 +123,4 @@ #endif + +