Index: trunk/MagicSoft/Mars/mtemp/meth/MAstroCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/MAstroCamera.cc (revision 5500) +++ (revision ) @@ -1,611 +1,0 @@ -/*====================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Thomas Bretz, 3/2004 -! Author(s): Sabrina Stark, 7/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MAstroCamera -// ============ -// -// A tools displaying stars from a catalog in the camera display. -// PRELIMINARY!! -// -// -// Usage: -// ------ -// For a usage example see macros/starfield.C -// -// To be able to reflect the star-light you need the geometry of the -// mirror and the distance of the plain screen. -// -// You can get the mirror geometry from a MC file and the distance of -// the screen from MGeomCam. -// -// -// Algorithm: -// ---------- -// The caluclation of the position of the reflection in the camera is -// done by: -// - Rotation of the star-field such that the camera is looking into -// the pointing direction -// - Calculation of the reflected star-light vector by calling -// MGeomMirror::GetReflection (which returns the point at which -// the vector hits the camera plain) -// - Depending on the draw-option you get each reflected point, the -// reflection on a virtual ideal mirror or the reflection on each -// individual mirror -// -// -// GUI: -// ---- -// * You can use the the cursor keys to change the pointing position -// and plus/minus to change the time by a quarter of an hour. -// -// ToDo: -// ----- -// * possibility to overwrite distance from mirror to screen -// * read the mirror geometry directly from the MC input file -// -///////////////////////////////////////////////////////////////////////////// -#include "MAstroCamera.h" - -#include // strerror -#include // ifstream - -#include // kKey_* - -#include // TH2D -#include // TMarker -#include // gPad - -#include "MLog.h" -#include "MLogManip.h" - -#include "MGeomCam.h" -#include "MGeomMirror.h" - -#include "MTime.h" -#include "MAstroSky2Local.h" -#include "MObservatory.h" - -#include "../mhist/MHCamera.h" // FIXME: This dependancy is very bad! - // HOW TO GET RID OF IT? Move MHCamera to mgeom? - -//#include "MStarLocalPos.h" - -ClassImp(MAstroCamera); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Create a virtual MGeomMirror which is in the center of the coordinate -// system and has a normal vector in z-direction. -// -MAstroCamera::MAstroCamera() : fGeom(0), fMirrors(0) -{ - fMirror0 = new MGeomMirror; - fMirror0->SetMirrorContent(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1); -} - -// -------------------------------------------------------------------------- -// -// Delete fGeom, fMirrors and the virtual 0-Mirror fMirror0 -// -MAstroCamera::~MAstroCamera() -{ - if (fGeom) - delete fGeom; - if (fMirrors) - delete fMirrors; - - delete fMirror0; - -} - -// -------------------------------------------------------------------------- -// -// Set a list of mirrors. The Mirrors must be of type MGeomMirror and -// stored in a TClonesArray -// -void MAstroCamera::SetMirrors(TClonesArray &arr) -{ - if (arr.GetClass()!=MGeomMirror::Class()) - { - cout << "ERROR - TClonesArray doesn't contain objects of type MGeomMirror... ignored." << endl; - return; - } - - const Int_t n = arr.GetSize(); - - if (!fMirrors) - fMirrors = new TClonesArray(MGeomMirror::Class(), n); - - fMirrors->ExpandCreate(n); - - for (int i=0; iExpandCreate(n); - continue; - } - - - Int_t id; - Float_t f, sx, sy, x, y, z, thetan, phin, xn, yn, zn; - - const Int_t n = sscanf(line.Data(), "%d %f %f %f %f %f %f %f %f %f %f %f", - &id, &f, &sx, &sy, &x, &y, &z, &thetan, - &phin, &xn, &yn, &zn); - if (n!=12) - continue; - - new ((*fMirrors)[id-1]) MGeomMirror; - ((MGeomMirror*)(*fMirrors)[id-1])->SetMirrorContent(id, f, sx, sy, x, y, z, thetan, phin, xn, yn, zn); - } -} - -// -------------------------------------------------------------------------- -// -// Set the camera geometry. The MGeomCam object is cloned. -// -void MAstroCamera::SetGeom(const MGeomCam &cam) -{ - if (fGeom) - delete fGeom; - - fGeom = (MGeomCam*)cam.Clone(); -} - -// -------------------------------------------------------------------------- -// -// Convert To Pad coordinates (see MAstroCatalog) -// -Int_t MAstroCamera::ConvertToPad(const TVector3 &w, TVector2 &v) const -{ - /* - --- Use this to plot the 'mean grid' instead of the grid of a - theoretical central mirror --- - - TVector3 spot; - const Int_t num = fConfig->GetNumMirror(); - for (int i=0; iGetMirror(i).GetReflection(w, fGeom->GetCameraDist())*1000; - spot *= 1./num; - */ - - const TVector3 spot = fMirror0->GetReflection(w, fGeom->GetCameraDist())*1000; - v.Set(spot(0), spot(1)); - - const Float_t max = fGeom->GetMaxRadius()*0.70; - return v.X()>-max && v.Y()>-max && v.X()GetListOfPrimitives()); - while ((o=Next())) - { - if (o->InheritsFrom(gROOT->GetClass(name))) - return o; - - if (o->InheritsFrom("TPad")) - if ((o = FindObjectInPad(name, (TVirtualPad*)o))) - return o; - } - return NULL; -} - -// -------------------------------------------------------------------------- -// -// Options: -// -// '*' Draw the mean of the reflections on all mirrors -// '.' Draw a dot for the reflection on each individual mirror -// 'h' To create a TH2D of the star-light which is displayed -// 'c' Use the underlaying MHCamera as histogram -// '0' Draw the reflection on a virtual perfect mirror -// '=' Draw '0' or '*' propotional to the star magnitude -// -// If the Pad contains an object MHCamera of type MHCamera it is used. -// Otherwise a new object is created. -// -void MAstroCamera::AddPrimitives(TString o) -{ - if (!fMirrors) - { - cout << "Missing Mirror data..." << endl; - } - if (!fObservatory) - { - cout << "Missing data (Observatory)..." << endl; - return; - } - if (!fTime) - { - cout << "Missing data (Time)..." << endl; - return; - } - - if (o.IsNull()) - o = "*."; - - const Bool_t hasnull = o.Contains("0", TString::kIgnoreCase); - const Bool_t hashist = o.Contains("h", TString::kIgnoreCase); - const Bool_t hasmean = o.Contains("*", TString::kIgnoreCase); - const Bool_t hasdot = o.Contains(".", TString::kIgnoreCase); - const Bool_t usecam = o.Contains("c", TString::kIgnoreCase); - const Bool_t resize = o.Contains("=", TString::kIgnoreCase); - - // Get camera - MHCamera *camera=(MHCamera*)FindObjectInPad("MHCamera", gPad); - if (camera) - { - if (!camera->GetGeometry() || camera->GetGeometry()->IsA()!=fGeom->IsA()) - camera->SetGeometry(*fGeom); - } - else - { - camera = new MHCamera(*fGeom); - camera->SetName("MHCamera"); - camera->SetStats(0); - camera->SetInvDeepBlueSeaPalette(); - camera->SetBit(kCanDelete); - camera->Draw(); - } - - camera->SetTitle(GetPadTitle()); - - gPad->cd(1); - - if (!usecam) - { - if (camera->GetEntries()==0) - camera->SetBit(MHCamera::kNoLegend); - } - else - { - camera->Reset(); - camera->SetYTitle("arb.cur"); - } - - TH2 *h=0; - if (hashist) - { - TH2F hist("","", 90, -650, 650, 90, -650, 650); - hist.SetMinimum(0); - h = (TH2*)hist.DrawCopy("samecont1"); - } - - const TRotation rot(GetGrid(kTRUE)); - - MVector3 *radec; - TIter Next(&fList); - - while ((radec=(MVector3*)Next())) - { - const Double_t mag = radec->Magnitude(); - if (mag>GetLimMag()) - continue; - - TVector3 star(*radec); - - // Rotate Star into telescope system - star *= rot; - - TVector3 mean; - - Int_t num = 0; - - MGeomMirror *mirror = 0; - TIter NextM(fMirrors); - while ((mirror=(MGeomMirror*)NextM())) - { - const TVector3 spot = mirror->GetReflection(star, fGeom->GetCameraDist())*1000; - - // calculate mean of all 'stars' hitting the camera plane - // by taking the sum of the intersection points between - // the light vector and the camera plane - mean += spot; - - if (hasdot) - { - TMarker *m=new TMarker(spot(0), spot(1), 1); - m->SetMarkerColor(kMagenta); - m->SetMarkerStyle(kDot); - fMapG.Add(m); - } - if (h) - h->Fill(spot(0), spot(1), pow(10, -mag/2.5)); - - if (usecam) - camera->Fill(spot(0), spot(1), pow(10, -mag/2.5)); - - num++; - } - - // transform meters into millimeters (camera display works with mm) - mean *= 1./num; - DrawStar(mean(0), mean(1), *radec, hasmean?kBlack:-1, Form("x=%.1fmm y=%.1fmm", mean(0), mean(1)), resize); - - if (hasnull) - { - TVector3 star(*radec); - star *= rot; - const TVector3 spot = fMirror0->GetReflection(star, fGeom->GetCameraDist())*1000; - DrawStar(spot(0), spot(1), *radec, hasmean?kBlack:-1, Form("x=%.1fmm y=%.1fmm", mean(0), mean(1)), resize); - // This can be used to get the abberation... - //cout << TMath::Hypot(spot(0), spot(1)) << " " << TMath::Hypot(mean(0)-spot(0), mean(1)-spot(1)) << endl; - } - } -} - -// -------------------------------------------------------------------------- -// -// Options: -// -// '*' Draw the mean of the reflections on all mirrors -// '.' Draw a dot for the reflection on each individual mirror -// 'h' To create a TH2D of the star-light which is displayed -// 'c' Use the underlaying MHCamera as histogram -// '0' Draw the reflection on a virtual perfect mirror -// -// If the Pad contains an object MHCamera of type MHCamera it is used. -// Otherwise a new object is created. -// - -void MAstroCamera::StarPosInCamera() -{ - if (!fMirrors) - { - cout << "Missing Mirror data..." << endl; - } - - const TRotation rot(GetGrid(kTRUE)); - - TIter Next(&fList); - MVector3 *v=0; - while ((v=(MVector3*)Next())) - { - const Double_t mag = v->Magnitude(); - if (mag>GetLimMag()) - continue; - - TVector3 star(*v); - - // Rotate Star into telescope system - star *= rot; - - TVector3 *mean = new TVector3; - - Int_t num = 0; - - // Aberration: - MGeomMirror *mirror = 0; - TIter NextM(fMirrors); - while ((mirror=(MGeomMirror*)NextM())) - { - const TVector3 spot = mirror->GetReflection(star, fGeom->GetCameraDist())*1000; - - *mean += spot; - num++; - } - - *mean *= 1./num; - fCatList.Add(mean); - } -} -/*void MAstroCamera::FillStarList(TList *list) -{ - list->SetOwner(); - list->Delete(); - - if (!fTime || !fObservatory || !fMirrors || !list) - { - cout << "Missing data..." << endl; - return; - } - - const MAstroSky2Local s2l(*fTime, *fObservatory); - const TRotation trans(AlignCoordinates(rot*fRaDec)); - - // Return the correct rotation matrix - const TRotation rot = trans*s2l; - - MVector3 *radec; - TIter Next(&fList); - - while ((radec=(MVector3*)Next())) - { - const Double_t mag = radec->Magnitude(); - - TVector3 mean; - TVector3 star(*radec); - star *= rot; - const TVector3 spot = fMirror0->GetReflection(star, fGeom->GetCameraDist())*1000; - - MStarLocalPos *starpos = new MStarLocalPos; - starpos->SetExpValues(mag,mean(0),mean(1)); - list->Add(starpos); - } - // For MAGIC the distance of the mean of the light distribution - // to the Mirror0 reflection of the star (Abberation) can be - // expressed as: dr = (0.0713 +/- 0.0002) * r = r/14.03 - // with r = hypot(mean(0), mean(1)) -} -*/ - -// ------------------------------------------------------------------------ -// -// Uses fRaDec as a reference point. -// -// Return dZd and dAz corresponding to the distance from x,y to fRaDec -// -// Before calling this function you should correct for abberation. In -// case of MAGIC you can do this by: -// x /= 1.0713; -// y /= 1.0713; -// -// x [mm]: x coordinate in the camera plane (assuming a perfect mirror) -// y [mm]: y coordinate in the camera plane (assuming a perfect mirror) -// -// We assume (0, 0) to be the center of the FOV -// -// dzd [deg]: Delta Zd -// daz [deg]: Delta Az -// -void MAstroCamera::GetDiffZdAz(Double_t x, Double_t y, Double_t &dzd, Double_t &daz) -{ - // Reflect the corrected pixel on a perfect mirror - TVector3 v(x, y, fGeom->GetCameraDist()*1000); - TVector3 spot = fMirror0->GetReflection(v); - - // Derotate into the local coordinate system - const MAstroSky2Local rot(*fTime, *fObservatory); - const TRotation align(AlignCoordinates(rot*fRaDec).Inverse()); - spot *= align; - - cout << "Zd="< -! Author(s): Javier López , 4/2004 -! Author(s): Sabrina Stark , 7/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MFindStars -// -///////////////////////////////////////////////////////////////////////////// -#include "MFindStars.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "MObservatory.h" -#include "MAstroCamera.h" -#include "MMcConfigRunHeader.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MHCamera.h" - -#include "MGeomCam.h" -#include "MGeomPix.h" -#include "MCameraDC.h" -#include "MTime.h" -#include "MReportDrive.h" -#include "MStarLocalCam.h" -#include "MStarLocalPos.h" - -#include "MParList.h" -#include "MTaskList.h" - -ClassImp(MFindStars); -using namespace std; - -const Float_t sqrt2 = sqrt(2.); -const Float_t sqrt3 = sqrt(3.); -const UInt_t lastInnerPixel = 396; - - -//______________________________________________________________________________ -// -// The 2D gaussian fucntion used to fit the spot of the star -// -static Double_t func(float x,float y,Double_t *par) -{ - Double_t value=par[0]+par[1]*exp(-(-(x-par[2])*sin(par[6]) - (y-par[4])*cos(par[6]))*(-(x-par[2])*sin(par[6]) - (y-par[4])*cos(par[6]))/(2*par[3]*par[3]))*exp(-((x-par[2])*cos(par[6]) - (y-par[4])*sin(par[6]))*((x-par[2])*cos(par[6]) - (y-par[4])*sin(par[6]))/(2*par[5]*par[5])); - - return value; -} - -//______________________________________________________________________________ -// -// Function used by Minuit to do the fit -// -static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) -{ - - 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=0; - - UInt_t usedPx=0; - for (UInt_t pixid=1; pixidlastInnerPixel?outerped:innerped); - errorz=(pixid>lastInnerPixel?outerrms:innerrms); - if (errorz > 0.0) - { - usedPx++; - delta = (z-func(x,y,par))/errorz; - chisq += delta*delta; - } - else - cerr << " TMinuit::fcn errorz[" << pixid << "] " << errorz << endl; - } - } - f = chisq; - - find->SetChisquare(chisq); - find->SetDegreesofFreedom(usedPx); -} - -MFindStars::MFindStars(const char *name, const char *title): - fGeomCam(NULL), fCurr(NULL), fTimeCurr(NULL), fDrive(NULL), fStars(NULL), fNumVar(7) -{ - fName = name ? name : "MFindStars"; - fTitle = title ? title : "Tool to find stars from DC Currents"; - - fNumIntegratedEvents=0; - fMaxNumIntegratedEvents = 10; - fRingInterest = 125.; //[mm] ~ 0.4 deg - fDCTailCut = 4; - - fPixelsUsed.Set(577); - fPixelsUsed.Reset((Char_t)kTRUE); - - //Fitting(Minuit) initialitation - const Float_t pixelSize = 31.5; //[mm] - fMinuitPrintOutLevel = -1; - - fVname = new TString[fNumVar]; - fVinit.Set(fNumVar); - fStep.Set(fNumVar); - fLimlo.Set(fNumVar); - fLimup.Set(fNumVar); - fFix.Set(fNumVar); - - fVname[0] = "background"; - fVinit[0] = fMaxNumIntegratedEvents; - fStep[0] = fVinit[0]/sqrt2; - fLimlo[0] = 0; - fLimup[0] = 4.*fMaxNumIntegratedEvents; - fFix[0] = 0; - - fVname[1] = "max"; - fVinit[1] = 10.*fMaxNumIntegratedEvents; - fStep[1] = fVinit[0]/sqrt2; - fLimlo[1] = fMinDCForStars; - fLimup[1] = 30.*fMaxNumIntegratedEvents; - fFix[1] = 0; - - fVname[2] = "meanx"; - fVinit[2] = 0.; - fStep[2] = fVinit[1]/sqrt2; - fLimlo[2] = -600.; - fLimup[2] = 600.; - fFix[2] = 0; - - fVname[3] = "sigmaminor"; - fVinit[3] = pixelSize; - fStep[3] = fVinit[2]/sqrt2; - fLimlo[3] = pixelSize/(2*sqrt3); - fLimup[3] = 500.; - fFix[3] = 0; - - fVname[4] = "meany"; - fVinit[4] = 0.; - fStep[4] = fVinit[3]/sqrt2; - fLimlo[4] = -600.; - fLimup[4] = 600.; - fFix[4] = 0; - - fVname[5] = "sigmamajor"; - fVinit[5] = pixelSize; - fStep[5] = fVinit[4]/sqrt2; - fLimlo[5] = pixelSize/(2*sqrt3); - fLimup[5] = 500.; - fFix[5] = 0; - - fVname[6] = "phi"; - fVinit[6] = 0.; - fStep[6] = fVinit[0]/sqrt2; - fLimlo[6] = 0.; - fLimup[6] = TMath::TwoPi(); - fFix[6] = 0; - - - fObjectFit = NULL; - // fMethod = "SIMPLEX"; - fMethod = "MIGRAD"; - // fMethod = "MINIMIZE"; - fNulloutput = kFALSE; - - // Set output level - // fLog->SetOutputLevel(3); // No dbg messages -} - -Int_t MFindStars::PreProcess(MParList *pList) -{ - - fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); - - if (!fGeomCam) - { - *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl; - return kFALSE; - } - - // Initialize camera display with the MGeomCam information - fDisplay.SetGeometry(*fGeomCam,"FindStarsDisplay","FindStarsDisplay"); - fDisplay.SetUsed(fPixelsUsed); - - fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC")); - - if (!fCurr) - { - *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl; - return kFALSE; - } - - fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents")); - - if (!fTimeCurr) - { - *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl; - return kFALSE; - } - - // fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive")); - - if (!fDrive) - { - *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl; - } - - - // Initalization of MAstroCamera - - fAstroCamera = new MAstroCamera(); - - // Current observatory - MObservatory magic1; - MMcConfigRunHeader *config=0; - MGeomCam *geom=0; - TString fname = "/dd/magLQE_3/gh/0/0/G_M0_00_0_550150_w0.root"; - TFile file(fname); - TTree *tree = (TTree*)file.Get("RunHeaders"); - tree->SetBranchAddress("MMcConfigRunHeader", &config); - if (tree->GetBranch("MGeomCam")) - tree->SetBranchAddress("MGeomCam", &geom); - tree->GetEntry(0); - - fAstroCamera->SetMirrors(*config->GetMirrors()); - fAstroCamera->SetGeom(*geom); - // Right Ascension [h] and declination [deg] of source - // Currently 'perfect' pointing is assumed - const Double_t ra = fAstro.Hms2Rad(fSrcRaHour,fSrcRaMin, fSrcRaSec); - const Double_t dec = fAstro.Dms2Rad(fSrcDecDeg, fSrcDecMin, fSrcDecSec); - cout << "----------------------------FindStars,pos = "<SetLimMag(6); - // Radius of FOV around the source position to load the stars - fAstroCamera->SetRadiusFOV(3); - // Source position - fAstroCamera->SetRaDec(ra, dec); - // Catalog to load (here: Bright Star Catalog V5) - fAstroCamera->ReadBSC("../catalogs/bsc5.dat"); - // Obersavatory and time to also get local coordinate information - fAstroCamera->SetObservatory(magic1); - // fRes contains the corrected source position - fRes = new TVector3; - - fStars = (MStarLocalCam*)pList->FindCreateObj(AddSerialNumber("MStarLocalCam")); - if (!fStars) - { - *fLog << err << AddSerialNumber("MStarLocalCam") << " cannot be created ... aborting" << endl; - return kFALSE; - } - - fMinDCForStars = 1.*fMaxNumIntegratedEvents; //[uA] - - // Initialize the TMinuit object - - TMinuit *gMinuit = new TMinuit(fNumVar); //initialize TMinuit with a maximum of params - gMinuit->SetFCN(fcn); - - Double_t arglist[2*fNumVar];//[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 MParList object pointer to allow mimuit to access internally - gMinuit->SetObjectFit(pList); - - return kTRUE; -} - -Int_t MFindStars::Process() -{ - UInt_t numPixels = fGeomCam->GetNumPixels(); - TArrayC origPixelsUsed; - origPixelsUsed.Set(numPixels); - if (fNumIntegratedEvents >= fMaxNumIntegratedEvents) - { - //Fist delete the previus stars in the list - fAstroCamera->GetCatList()->Delete(); - // Time for which to get the star position - fAstroCamera->SetTime(*fTimeCurr); // FIX ME: time should be averaged - //Clone the catalog due to the validity range of the instance - TObject *o = fAstroCamera->Clone(); - o->SetBit(kCanDelete); - o->Draw("mirror"); - - fAstroCamera->StarPosInCamera(); - - if (fDrive) - { -// Float_t ra = fDrive->GetRa(); -// Float_t dec = fDrive->GetDec(); - -// fAstro.SetRaDec(ra, dec); -// fAstro.SetGuiActive(); - -// fAstro.FillStarList(); - } - else - { - //Fist delete the previus stars in the list - fStars->GetList()->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); - } - - //Get the corrected source position - TVector3 *res = new TVector3; - if (CorrSourcePos(res)) - fRes = res; - else - fRes->SetXYZ(-1000, -1000, -1000); - //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=0) *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <GetNumPixels(); - 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); - - for (UInt_t pix=1; pix<=lastInnerPixel; pix++) - dchist[0]->Fill(fDisplay.GetBinContent(pix+1)); - for (UInt_t pix=lastInnerPixel+1; pixFill(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; -} - -Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix) -{ - UInt_t numPixels = fGeomCam->GetNumPixels(); - -// Find the two close pixels with the maximun dc - UInt_t maxPixIdx[2]; - - maxDC = 0; - - for (UInt_t pix=1; pix maxDC*2) - { - if(dc[0]>=dc[1]) - { - maxPixIdx[0] = pix; - maxPixIdx[1] = swneighbor; - maxDC = dc[0]; - } - else - { - maxPixIdx[1] = pix; - maxPixIdx[0] = swneighbor; - maxDC = dc[1]; - } - } - } - } - } - } - - if (maxDC == 0) - { - *fLog << warn << " No found pixels with maximum dc" << endl; - return kFALSE; - } - - maxPix = (*fGeomCam)[maxPixIdx[0]]; - - 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; -} - -Bool_t MFindStars::FindStar(MStarLocalPos* star) -{ - - UInt_t numPixels = fGeomCam->GetNumPixels(); - Float_t innerped = fStars->GetInnerPedestalDC(); - Float_t outerped = fStars->GetOuterPedestalDC(); - - TArrayC origPixelsUsed; - origPixelsUsed.Set(numPixels); - - for (UInt_t pix=1; pixGetXExp(); - Float_t expY = star->GetYExp(); - - Float_t max=0; - UInt_t pixmax=0; - Float_t meanX=0; - Float_t meanY=0; - Float_t meanSqX=0; - Float_t meanSqY=0; - Float_t sumCharge=0; - UInt_t usedInnerPx=0; - UInt_t usedOuterPx=0; - - const Float_t meanPresition = 3.; //[mm] - const UInt_t maxNumIterations = 10; - UInt_t numIterations = 0; - - do - { - // First define a area of interest around the expected position of the star - 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; - 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; - - -// fit the star spot using TMinuit - - - for (UInt_t pix=1; pix=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] = 0; - fVinit[1] = max; - fLimlo[1] = fMinDCForStars-(pixmax>lastInnerPixel?outerped:innerped); - fLimup[1] = 30*fMaxNumIntegratedEvents-(pixmax>lastInnerPixel?outerped:innerped); - fVinit[2] = meanX; - fVinit[3] = rmsX; - fVinit[4] = meanY; - fVinit[5] = rmsY; - fVinit[6] = 0; - //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; - } - // - - // ------------------------------------------- - // 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(); - -// Now ready for minimization step - arglist[0] = 500; - arglist[1] = 1.; - gMinuit->mnexcm(fMethod, arglist ,2,ierflg); - - clock.Stop(); - - if(fMinuitPrintOutLevel>=0) - { - if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;; - clock.Print(); - } - - Double_t integratedCharge; - Double_t bkFit, bkFitError; - Double_t maxFit, maxFitError; - Double_t meanXFit, meanXFitError; - Double_t sigmaMinorAxis, sigmaMinorAxisError; - Double_t meanYFit, meanYFitError; - Double_t sigmaMajorAxis, sigmaMajorAxisError; - Double_t phiFit, phiFitError; - Float_t chisquare = GetChisquare(); - Int_t dregrees = GetDegreesofFreedom()-fNumVar; - - if (!ierflg) - { - gMinuit->GetParameter(0,bkFit, bkFitError); - gMinuit->GetParameter(1,maxFit, maxFitError); - gMinuit->GetParameter(2,meanXFit,meanXFitError); - gMinuit->GetParameter(3,sigmaMinorAxis,sigmaMinorAxisError); - gMinuit->GetParameter(4,meanYFit,meanYFitError); - gMinuit->GetParameter(5,sigmaMajorAxis,sigmaMajorAxisError); - gMinuit->GetParameter(6,phiFit,phiFitError); - - //FIXME: Do the integral properlly - integratedCharge = 0.; - - - } - else - { - bkFit = 0.; - maxFit = 0.; - meanXFit = 0.; - sigmaMinorAxis = 0.; - meanYFit = 0.; - sigmaMajorAxis = 0.; - integratedCharge = 0.; - phiFit = 0.; - - *fLog << err << "TMinuit::Call error " << ierflg << endl; - } - - star->SetFitValues(integratedCharge,bkFit,maxFit,meanXFit,meanYFit,sigmaMinorAxis,sigmaMajorAxis,phiFit,chisquare,dregrees); - - // reset the display to the starting values - fDisplay.SetUsed(origPixelsUsed); - - if (ierflg) - return kCONTINUE; - return kTRUE; -} - -Bool_t MFindStars::ShadowStar(MStarLocalPos* star) -{ - UInt_t numPixels = fGeomCam->GetNumPixels(); - -// Define an area around the star which will be set unused. - UInt_t shadowPx=0; - for (UInt_t pix=1; pixGetMeanX(); - Float_t starYpos = star->GetMeanY(); - - Float_t starSize = 3*star->GetSigmaMajorAxis(); - - Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+ - (pixYpos-starYpos)*(pixYpos-starYpos)); - - if (dist > starSize && fDisplay.IsUsed(pix)) - fPixelsUsed[pix]=(Char_t)kTRUE; - else - { - fPixelsUsed[pix]=(Char_t)kFALSE; - shadowPx++; - } - } - - if (fMinuitPrintOutLevel>=0) *fLog << dbg << " shadowPx " << shadowPx << endl; - - fDisplay.SetUsed(fPixelsUsed); - - return kTRUE; -} - -Bool_t MFindStars::CorrSourcePos(TVector3* srcpos) -{ - const Int_t numstar = fStars->GetList()->GetSize(); - if (numstar < 2) - { - cout <<" Less than 2 stars in list. Triangulation not possible."<GetList()); - MStarLocalPos* star; - while ((star=(MStarLocalPos*)Next())) - { - starmag[i] = star->GetMagCalc(); - starposx[i] = star->GetMeanX(); - starposy[i] = star->GetMeanY(); - i +=1; - } - for ( i = 0; i < numstar; i++) - { - starposcatx[i] = -1000; - starposcaty[i] = -1000; - for( Int_t dist = 10; dist < maxdist; dist += 10) - { - TIter Next1(fAstroCamera->GetCatList()); - TVector3 *starcat; - while ((starcat=(TVector3*)Next1()) && dist < maxdist) - { - if ((starcat->X()-starposx[i])*(starcat->X()-starposx[i])+(starcat->Y()-starposy[i])*(starcat->Y()-starposy[i]) < dist*dist) - { - starposcatx[i] = starcat->X(); - starposcaty[i] = starcat->Y(); - dist = maxdist; - matches ++; - } - } - } - } - // take the 2 'best' determined stars: Assumpion that the bright stars are better. - Float_t maxmag1 = 0, maxmag2 = 0; - Float_t minchi1 = 0, minchi2 = 0; - Int_t magi1 = -1, magi2 = -1; - Int_t chii1 = -1, chii2 = -1; - for ( i = 0; i < numstar; i++) - { - if (starposcatx[i] != -1000 && starposcaty[i] != -1000) - { - if (starmag[i] > maxmag2) - { - maxmag2 = starmag[i]; - magi2 = i; - } - if (maxmag2 > maxmag1) - { - maxmag2 = maxmag1; - maxmag1 = starmag[i]; - magi2 = magi1; - magi1 = i; - } - } - } - - if (matches < 2) - { - cout << " Could not find 2 stars in the camera matching with the stars in the catalog. Triangulation not possible."<Set(0.,0.); // Assume correct pointing. - star1->Set(starposcatx[magi1],starposcaty[magi1]); - star2->Set(starposcatx[magi2],starposcaty[magi2]); - - DistBetweenStars(star1,star2,source,dist); - - Float_t dx = starposx[magi1] - starposx[magi2]; - Float_t delta = acos(dx/dist->Y()); - Float_t alpha = acos((dist->Y()*dist->Y()+dist->Z()*dist->Z()-dist->X()*dist->X())/2/dist->Y()/dist->Y()); - Float_t sigma = alpha - delta; - - if (source->X()-starposx[magi1] > 0) - srcpos->SetX(starposx[magi1] + dist->Z()*abs(cos(sigma))); - else - srcpos->SetX(starposx[magi1] - dist->Z()*abs(cos(sigma))); - if (source->Y()-starposy[magi1] > 0) - srcpos->SetY(starposy[magi1] + dist->Z()*abs(sin(sigma))); - else - srcpos->SetY(starposy[magi1] - dist->Z()*abs(sin(sigma))); - return kTRUE; -} - - -Bool_t MFindStars::DistBetweenStars(TVector2 *star1, TVector2 *star2, TVector2 *source, TVector3 *dist) -{ - dist->SetX(sqrt(abs((source->X()-star2->X())*(source->X()-star2->X())+(source->Y()-star2->Y())*(source->Y()-star2->Y())))); - dist->SetY(sqrt(abs((star1->X()-star2->X())*(star1->X()-star2->X())+(star1->Y()-star2->Y())*(star1->Y()-star2->Y())))); - dist->SetZ(sqrt(abs((source->X()-star1->X())*(source->X()-star1->X())+(source->Y()-star1->Y())*(source->Y()-star1->Y())))); - return kTRUE; - -} Index: trunk/MagicSoft/Mars/mtemp/meth/MFindStars.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/MFindStars.h (revision 5500) +++ (revision ) @@ -1,134 +1,0 @@ -#ifndef MARS_MFindStars -#define MARS_MFindStars - -#ifndef ROOT_TArrayS -#include -#endif - -#ifndef ROOT_TArrayC -#include -#endif - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -#ifndef MARS_MHCamera -#include "MHCamera.h" -#endif - -#ifndef MARS_MAstroCamera -#include "MAstroCamera.h" -#endif - -#ifndef MARS_MAstro -#include "MAstro.h" -#endif - -class MGeomCam; -class MGeomPix; -class MCameraDC; -class MTime; -class MReportDrive; -class MStarLocalCam; -class MStarLocalPos; - -class MFindStars : public MTask -{ - -private: - - MGeomCam *fGeomCam; - MCameraDC *fCurr; - MTime *fTimeCurr; - MReportDrive *fDrive; - MStarLocalCam *fStars; - - MAstroCamera *fAstroCamera; - MAstro fAstro; - TArrayC fPixelsUsed; - MHCamera fDisplay; - - UInt_t fMaxNumIntegratedEvents; - UInt_t fNumIntegratedEvents; - - Float_t fRingInterest; //[mm] - Float_t fMinDCForStars; //[uA] - Float_t fDCTailCut; //[uA] - - // Source position: - Int_t fSrcRaHour; - Int_t fSrcRaMin; - Float_t fSrcRaSec; - - Int_t fSrcDecDeg; - Int_t fSrcDecMin; - Float_t fSrcDecSec; - - TVector3 *fRes; - - //Fitting(Minuit) variables - const Int_t fNumVar; - Float_t fTempChisquare; - Int_t fTempDegreesofFreedom; - Int_t fMinuitPrintOutLevel; - - TString *fVname; - TArrayD fVinit; - TArrayD fStep; - TArrayD fLimlo; - TArrayD fLimup; - TArrayI fFix; - TObject *fObjectFit; - TString fMethod; - Bool_t fNulloutput; - - Bool_t DCPedestalCalc(); - Bool_t FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix); - Bool_t FindStar(MStarLocalPos* star); - Bool_t ShadowStar(MStarLocalPos* star); - Bool_t CorrSourcePos(TVector3* srcpos); - Bool_t DistBetweenStars(TVector2* star1, TVector2* star2, TVector2* source,TVector3* dist); - - public: - - MFindStars(const char *name=NULL, const char *title=NULL); - - Int_t PreProcess(MParList *pList); - 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;} - - void SetSourceRaH(Int_t hour) {fSrcRaHour=hour;} - void SetSourceRaM(Int_t min) {fSrcRaMin=min;} - void SetSourceRaS(Float_t sec) {fSrcRaSec=sec;} - - void SetSourceDecD(Int_t deg) {fSrcDecDeg=deg;} - void SetSourceDecM(Int_t min) {fSrcDecMin=min;} - void SetSourceDecS(Float_t sec) {fSrcDecSec=sec;} - - - //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;} - TList *GetCatalogList() { return fAstroCamera->GetCatList(); } - TVector3 *GetCorrSrcPos() {return fRes;} - - - ClassDef(MFindStars, 0) // Tool to find stars from DC Currents -}; - -#endif Index: trunk/MagicSoft/Mars/mtemp/meth/MStarLocalPos.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/MStarLocalPos.cc (revision 5500) +++ (revision ) @@ -1,196 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Javier López , 4/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ -#include "MStarLocalPos.h" - -#include - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MStarLocalPos); - -using namespace std; - -MStarLocalPos::MStarLocalPos(const char *name, const char *title) -{ - - fName = name ? name : "MStarLocalPos"; - fTitle = title ? title : ""; - - Reset(); -} - -void MStarLocalPos::Reset() -{ - - //Expected position on camera - fMagExp = 0.; - fXExp = 0.; - fYExp = 0.; - - //Info from calculation - - fMagCalc = 0.; - fMaxCalc = 0.; - fMeanXCalc = 0.; - fMeanYCalc = 0.; - fSigmaMinorAxisCalc = 0.; - fSigmaMajorAxisCalc = 0.; - - //Info from fit - - fMagFit = 0.; - fBkFit = 0.; - fMaxFit = 0.; - fMeanXFit = 0.; - fMeanYFit = 0.; - fSigmaMinorAxisFit = 0.; - fSigmaMajorAxisFit = 0.; - fPhiFit = 0.; - fChiSquare = 0.; - fNdof = 1; - -} - -void MStarLocalPos::SetExpValues(Float_t mag, Float_t x, Float_t y) -{ - fMagExp = mag; - fXExp = x; - fYExp = y; -} - -void MStarLocalPos::SetCalcValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis) -{ - fMagCalc = mag; - fMaxCalc = max; - fMeanXCalc = x; - fMeanYCalc = y; - fSigmaMinorAxisCalc = sigmaMinorAxis; - fSigmaMajorAxisCalc = sigmaMajorAxis; -} - -void MStarLocalPos::SetFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, Float_t chiSquare, Int_t ndof) -{ - fMagFit = mag; - fMaxFit = max; - fMeanXFit = x; - fMeanYFit = y; - fSigmaMinorAxisFit = sigmaMinorAxis; - fSigmaMajorAxisFit = sigmaMajorAxis; - fChiSquare = chiSquare; - fNdof = ndof; -} - -void MStarLocalPos::SetFitValues(Float_t mag, Float_t bk, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, Float_t phi, Float_t chiSquare, Int_t ndof) -{ - fMagFit = mag; - fBkFit = bk; - fMaxFit = max; - fMeanXFit = x; - fMeanYFit = y; - fSigmaMinorAxisFit = sigmaMinorAxis; - fSigmaMajorAxisFit = sigmaMajorAxis; - fPhiFit = phi; - fChiSquare = chiSquare; - fNdof = ndof; -} - -// -------------------------------------------------------------------------- -// -// Paint the ellipse corresponding to the parameters -// -void MStarLocalPos::Paint(Option_t *opt) -{ - //Print a cross in the expected position - - if (fSigmaMinorAxisCalc>0. && fSigmaMajorAxisCalc>0.) - { - TEllipse ecalc(fMeanXCalc, fMeanYCalc, fSigmaMinorAxisCalc, fSigmaMajorAxisCalc, 0, 360, 0); - ecalc.SetLineWidth(2); - ecalc.SetLineColor(kRed); - ecalc.Paint(); - } - - if (fSigmaMinorAxisFit>0. || fSigmaMajorAxisFit>0.) - { - TEllipse efit(fMeanXFit, fMeanYFit, fSigmaMinorAxisFit, fSigmaMajorAxisFit, 0, 360, fPhiFit); - efit.SetLineWidth(2); - efit.SetLineColor(kBlack); - efit.Paint(); - } -} - -void MStarLocalPos::Print(Option_t *opt) const -{ - TString o = opt; - - if (o.Contains("mag", TString::kIgnoreCase) || opt == NULL) - { - *fLog << inf << "Star maginitude:" << endl; - *fLog << inf << " Expected \t" << setw(4) << fMagExp << endl; - *fLog << inf << " Calcultated \t " << setw(4) << fMagCalc << endl; - *fLog << inf << " Fitted \t " << setw(4) << fMagFit << endl; - } - - if (o.Contains("bk", TString::kIgnoreCase) || opt == NULL) - { - *fLog << inf << "Background:" << endl; - *fLog << inf << " Fitted \t " << setw(4) << fBkFit << " uA" << endl; - } - - if (o.Contains("max", TString::kIgnoreCase) || opt == NULL) - { - *fLog << inf << "Star Maximum:" << endl; - *fLog << inf << " Calcultated \t " << setw(4) << fMaxCalc << " uA" << endl; - *fLog << inf << " Fitted \t " << setw(4) << fMaxFit << " uA" << endl; - } - - if (o.Contains("pos", TString::kIgnoreCase) || opt == NULL) - { - *fLog << inf << "Star position:" << endl; - *fLog << inf << " Expected \t X " << setw(4) << fXExp << " mm \tY " << setw(4) << fYExp << " mm" << endl; - *fLog << inf << " Calcultated \t X " << setw(4) << fMeanXCalc << " mm \tY " << setw(4) << fMeanYCalc << " mm" << endl; - *fLog << inf << " Fitted \t X " << setw(4) << fMeanXFit << " mm \tY " << setw(4) << fMeanYFit << " mm" << endl; - } - - if (o.Contains("rot", TString::kIgnoreCase) || opt == NULL) - { - *fLog << inf << "Rotation of ellipse:" << endl; - *fLog << inf << " Fitted \t " << setw(4) << fPhiFit << " rad " << endl; - } - - if (o.Contains("siz", TString::kIgnoreCase) || opt == NULL) - { - *fLog << inf << "Star size:" << endl; - *fLog << inf << " Calcultated \t X " << setw(4) << fSigmaMinorAxisCalc << " mm \tY " << setw(4) << fSigmaMajorAxisCalc << " mm" << endl; - *fLog << inf << " Fitted \t X " << setw(4) << fSigmaMinorAxisFit << " mm \tY " << setw(4) << fSigmaMajorAxisFit << " mm" << endl; - } - - if (o.Contains("chi", TString::kIgnoreCase) || opt == NULL) - { - *fLog << inf << "Star Fit Quality:" << endl; - *fLog << inf << " ChiSquare/Ndof \t " << setw(3) << fChiSquare << "/" << fNdof << endl; - } -} - Index: trunk/MagicSoft/Mars/mtemp/meth/MStarLocalPos.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/MStarLocalPos.h (revision 5500) +++ (revision ) @@ -1,86 +1,0 @@ -#ifndef MARS_MStarLocalPos -#define MARS_MStarLocalPos - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -class MStarLocalPos : public MParContainer -{ -private: - - //Expected position on camera - - Float_t fMagExp; - Float_t fXExp; //[mm] - Float_t fYExp; //[mm] - - //Info from calculation - - Float_t fMagCalc; - Float_t fMaxCalc; //[uA] - Float_t fMeanXCalc; //[mm] - Float_t fMeanYCalc; //[mm] - Float_t fSigmaMinorAxisCalc; //[mm] - Float_t fSigmaMajorAxisCalc; //[mm] - - //Info from fit - - Float_t fMagFit; - Float_t fBkFit; //[uA] - Float_t fMaxFit; //[uA] - Float_t fMeanXFit; //[mm] - Float_t fMeanYFit; //[mm] - Float_t fSigmaMinorAxisFit; //[mm] - Float_t fSigmaMajorAxisFit; //[mm] - Float_t fPhiFit; //[rad] - Float_t fChiSquare; - Int_t fNdof; - -public: - - MStarLocalPos(const char *name=NULL, const char *title=NULL); - //~MStarLocalPos(); - - Float_t GetMagExp() {return fMagExp;} - Float_t GetXExp() {return fXExp;} - Float_t GetYExp() {return fYExp;} - - Float_t GetMagCalc() {return fMagCalc;} - Float_t GetMaxCalc() {return fMaxCalc;} - Float_t GetMeanXCalc() {return fMeanXCalc;} - Float_t GetMeanYCalc() {return fMeanYCalc;} - Float_t GetSigmaMinorAxisCalc() {return fSigmaMinorAxisCalc;} - Float_t GetSigmaMajorAxisCalc() {return fSigmaMajorAxisCalc;} - - Float_t GetMagFit() {return fMagFit;} - Float_t GetBkFit() {return fBkFit;} - Float_t GetMaxFit() {return fMaxFit;} - Float_t GetMeanXFit() {return fMeanXFit;} - Float_t GetMeanYFit() {return fMeanYFit;} - Float_t GetSigmaMinorAxisFit() {return fSigmaMinorAxisFit;} - Float_t GetSigmaMajorAxisFit() {return fSigmaMajorAxisFit;} - Float_t GetPhiFit() {return fPhiFit;} - Float_t GetChiSquare() {return fChiSquare;} - UInt_t GetNdof() {return fNdof;} - Float_t GetChiSquareNdof() {return fChiSquare/fNdof;} - - Float_t GetMeanX() {return fMeanXFit!=0?fMeanXFit:fMeanXCalc;} - Float_t GetMeanY() {return fMeanXFit!=0?fMeanYFit:fMeanYCalc;} - Float_t GetSigmaMinorAxis() {return fSigmaMinorAxisFit!=0?fSigmaMinorAxisFit:fSigmaMinorAxisCalc;} - Float_t GetSigmaMajorAxis() {return fSigmaMajorAxisFit!=0?fSigmaMajorAxisFit:fSigmaMajorAxisCalc;} - - void Reset(); - - void SetExpValues(Float_t mag, Float_t x, Float_t y); - void SetCalcValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis); - void SetFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, Float_t chi, Int_t ndof); - void SetFitValues(Float_t mag, Float_t bk, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, Float_t phi, Float_t chi, Int_t ndof); - - void Paint(Option_t *opt=NULL); - void Print(Option_t *opt=NULL) const; - - ClassDef(MStarLocalPos, 1) // Container that holds the star information in the PMT camera -}; - -#endif