Index: trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/Changelog (revision 4413) +++ trunk/MagicSoft/Mars/mtemp/Changelog (revision 4415) @@ -19,5 +19,34 @@ -*-*- END OF LINE -*-*- - 2004/07/20: David Paneque + 2004/07/22: Sabrina Stark + + * /mtemp/meth/MFindStars.[h,cc] + - new fit function, incl. constant background and rotation of + 2d gauss + - initialization of MAstroCamera + - new functions: + - CorrSourcePos -> calculates a corrected + source position using the images of two stars + - DistBetweenStars -> calculates the distances between + two stars and the source + + * /mtemp/meth/MAstroCamera.[h,cc] + - TList *fCatList: list with star positions from catalog + - new functions: + - StarPosInCamera (similar to FillStarList, but + incl. aberration) + - GetCatList -> Access function to fCatList + + * /mtemp/meth/MStarLocalPos.[h,cc] + - in function SetFitValue -> background included + + + * /mtemp/meth/CorrectSrcPos.C + - macro to calculate a corrected source position using + the images of two stars, these positions are written + to an output file + + + 2004/07/20: David Paneque * /mtemp/mmpi/SupercutsONOFFClasses directory containing the Index: trunk/MagicSoft/Mars/mtemp/meth/CorrectSrcPos.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/CorrectSrcPos.C (revision 4415) +++ trunk/MagicSoft/Mars/mtemp/meth/CorrectSrcPos.C (revision 4415) @@ -0,0 +1,359 @@ +/* ======================================================================== *\ +! +! * +! * 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): Luisa Sabrina Stark Schneebeli, 07/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! This macro is partially derived from findstars.C +\* ======================================================================== */ + + +Bool_t HandleInput() +{ + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + while (1) + { + // + // While reading the input process gui events asynchronously + // + timer.TurnOn(); + TString input = Getline("Type 'q' to exit, to go on: "); + timer.TurnOff(); + + if (input=="q\n") + return kFALSE; + + if (input=="\n") + return kTRUE; + }; + + return kFALSE; +} + + +void CorrectSrcPos(const TString filename="dc_2004_04_21_22_*_Mrk421.root", const TString directory="/scratch/Period016/cacodata/2004_04_22/", const UInt_t numEvents = 0) +{ + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + + MGeomCamMagic geomcam; + MCameraDC dccam; + MStarLocalCam starcam; + + plist.AddToList(&geomcam); + plist.AddToList(&dccam); + plist.AddToList(&starcam); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + + // Reads the trees of the root file and the analysed branches + MReadReports read; + read.AddTree("Currents"); + read.AddFile(directory+filename); // after the reading of the trees!!! + read.AddToBranchList("MReportCurrents.*"); + + MGeomApply geomapl; + TString continuoslightfile = + "/scratch/Period016/cacodata/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; + MCalibrateDC dccal; + dccal.SetFileName(continuoslightfile); + + + const Int_t numblind = 12; + const Short_t x[numblind] = { 8, 27, + 507, 508, 509, 510, 511, + 543, 559, 560, 561, 567}; + const TArrayS blindpixels(numblind,(Short_t*)x); + Float_t ringinterest = 100; //[mm] + Float_t tailcut = 3.5; + UInt_t integratedevents = 5; + + MFindStars findstars; + // findstars.SetBlindPixels(blindpixels); + findstars.SetRingInterest(ringinterest); + findstars.SetDCTailCut(tailcut); + findstars.SetNumIntegratedEvents(integratedevents); + findstars.SetMinuitPrintOutLevel(0); + + findstars.SetSourceRaH(11); // Mrk421 + findstars.SetSourceRaM(4); + findstars.SetSourceRaS(26.0); + findstars.SetSourceDecD(38); + findstars.SetSourceDecM(12); + findstars.SetSourceDecS(36.0); + + + // prints + MPrint pdc("MCameraDC"); + MPrint pstar("MStarLocalCam"); + + tlist.AddToList(&geomapl); + tlist.AddToList(&read); + tlist.AddToList(&dccal); + // tlist.AddToList(&pdc, "Currents"); + tlist.AddToList(&findstars, "Currents"); + // tlist.AddToList(&pstar, "Currents"); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + const Int_t numpar = 3; + Float_t starpos[10*numpar]; + Float_t starposcat[10*numpar]; + + if (numEvents > 0) + { + if (!evtloop.Eventloop(numEvents)) + return; + } + else + { + if (!evtloop.PreProcess()) + return; + + MHCamera display(geomcam); + display.SetPrettyPalette(); + + TCanvas *c = new TCanvas("c", "FindStars", 100, 100,600,600); + + TH2D *h0 = new TH2D("h0","",1,-2,2,1,-2,2); + h0->GetXaxis()->SetTitle("(deg)"); + h0->GetYaxis()->SetTitle("(deg)"); + h0->Draw(); + + display.Draw(); + gPad->cd(1); + //display.Draw("same"); + starcam.Draw(); + + UShort_t year; + Byte_t month, day; + + MTime *evttime = (MTime*)plist->FindObject("MTimeCurrents"); + + TList *starcatalog = new TList; + TList *stardc = new TList; + TVector3 *res = new TVector3; + + Int_t maxdist = 100; + UInt_t numevents=0; + Int_t numEvt=0; + const Int_t maxnumEvt = 10000; + Int_t col = 2; + const Double_t kDeg2mm = 189/0.6; + + Double_t resx[maxnumEvt]; + Double_t resy[maxnumEvt]; + Double_t time[maxnumEvt]; + + +// +// Find in the file name the date, run and project name + + Size_t pos = filename.Last('/'); + TString iRun = TString(filename(pos+24,5)); + TString iYear = TString(filename(pos+4,4)); + TString iMonth = TString(filename(pos+9,2)); + TString iDay = TString(filename(pos+12,2)); + + TString iHour = TString(filename(pos+15,2)); + TString iMin = TString(filename(pos+18,2)); + TString iSec = TString(filename(pos+21,2)); + + Size_t poslast = filename.Last('.'); + Size_t posfirst = poslast-1; + while (filename[posfirst] != '_') + posfirst--; + + TString iSource = TString(filename(posfirst+1,poslast-posfirst-1)); + char str[100]; + char str2[100]; + sprintf(str,"Date %s %s %s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); + sprintf(str2,"CorrectSrcPos_%s%s%s_%s_%s.txt",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); + + UInt_t evtyear, evtmonth, evtday, evthour, evtmin, evtsec; + UShort_t evtmsec; + + ofstream fposout1; + // fposout1.open (str2, ofstream::app); + fposout1.open (str2, ofstream::trunc); + + while (tlist.Process()) + { + if (numevents >= integratedevents) + { + evttime->GetTime(evthour, evtmin, evtsec,evtmsec); + evttime->Print(); + evthour = evttime->Hour(); + evtmin = evttime->Min(); + evtsec = evttime->Sec(); + + // display.SetCamContent(dccam); + display.SetCamContent(findstars.GetDisplay()); + gPad->Modified(); + gPad->Update(); + Int_t numStars = starcam.GetNumStars(); + starcam.Print("bkmaxposrotsizechi"); + + res = findstars.GetCorrSrcPos(); + cout<<"Corrected source position = "<X()<<" , "<Y()<X(); + resy[numEvt] = res->Y(); + time[numEvt] = numEvt; + numEvt ++; + + // write the corrected source position and the star positions to file + if (!fposout1) + cout << "Could not open output file to write the center of the fit." << endl; + else{ + for ( Int_t i = 0; i < numStars; i++) + { + if (i == 0){ + fposout1 << str<<" "<X()<< " " << res->Y(); + } + fposout1 <<" "<< starcam[i].GetMeanX()<< " " << starcam[i].GetMeanY(); + } + } + fposout1 << endl; + numevents = 0; + + // Remove the comments if you want to go through the file + // event-by-event: + //if (!HandleInput()) + // break; + } + numevents++; + } + fposout1.close(); + + // Print source position + gROOT->Reset(); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetOptStat(00000000); + + TCanvas *c1 = new TCanvas("c1", "SourcePosition", 200, 200,600,1200); + c1->Divide(1,3); + + TMath math; + + Double_t minresx, maxresx; + minresx = resx[math.LocMin(numEvt,resx)]; + maxresx = resx[math.LocMax(numEvt,resx)]; + + Double_t minresy, maxresy; + minresy = resy[math.LocMin(numEvt,resy)]; + maxresy = resy[math.LocMax(numEvt,resy)]; + + Double_t minres, maxres; + minres = math.Min(minresx,minresy); + maxres = math.Max(maxresx,maxresy); + + Double_t diff; + diff = maxres - minres; + diff = 0.1*diff; + minres = minres - diff; + maxres = maxres + diff; + + Double_t mintime, maxtime; + mintime = time[math.LocMin(numEvt,time)]; + maxtime = time[math.LocMax(numEvt,time)]; + + c1->cd(1); + TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minres,maxres); + h1->GetXaxis()->SetTitle("(evt.nr)"); + h1->GetYaxis()->SetTitle("(mm)"); + h1->Draw(); + + TGraph *grtimeevolx = new TGraph(numEvt,time,resx); + grtimeevolx->SetMarkerColor(215); + grtimeevolx->SetMarkerStyle(20); + grtimeevolx->SetMarkerSize (0.4); + grtimeevolx->Draw("P"); + + TGraph *grtimeevoly = new TGraph(numEvt,time,resy); + grtimeevoly->SetMarkerColor(6); + grtimeevoly->SetMarkerStyle(24); + grtimeevoly->SetMarkerSize (0.4); + grtimeevoly->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevolx,"src pos x","P"); + legxy.AddEntry(grtimeevoly,"src pos y","P"); + legxy.Draw(); + + c1->cd(2); + TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minresx,maxresx); + h2->GetXaxis()->SetTitle("(evt.nr)"); + h2->GetYaxis()->SetTitle("(mm)"); + h2->Draw(); + + TGraph *grtimeevolx = new TGraph(numEvt,time,resx); + grtimeevolx->SetMarkerColor(215); + grtimeevolx->SetMarkerStyle(20); + grtimeevolx->SetMarkerSize (0.4); + grtimeevolx->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevolx,"src pos x","P"); + legxy.Draw(); + + c1->cd(3); + TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minresy,maxresy); + h3->GetXaxis()->SetTitle("(evt.nr)"); + h3->GetYaxis()->SetTitle("(mm)"); + h3->Draw(); + + TGraph *grtimeevoly = new TGraph(numEvt,time,resy); + grtimeevoly->SetMarkerColor(6); + grtimeevoly->SetMarkerStyle(24); + grtimeevoly->SetMarkerSize (0.4); + grtimeevoly->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevoly,"src pos y","P"); + legxy.Draw(); + + evtloop.PostProcess(); + } + + tlist.PrintStatistics(); + +} + + Index: trunk/MagicSoft/Mars/mtemp/meth/MAstroCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/MAstroCamera.cc (revision 4415) +++ trunk/MagicSoft/Mars/mtemp/meth/MAstroCamera.cc (revision 4415) @@ -0,0 +1,611 @@ +/*====================================================================== *\ +! +! * +! * 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); + // Magnitude up to which the stars are loaded from the catalog + fAstroCamera->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; +} + + +void 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()))); + +} Index: trunk/MagicSoft/Mars/mtemp/meth/MFindStars.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/MFindStars.h (revision 4415) +++ trunk/MagicSoft/Mars/mtemp/meth/MFindStars.h (revision 4415) @@ -0,0 +1,134 @@ +#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); + void 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 4415) +++ trunk/MagicSoft/Mars/mtemp/meth/MStarLocalPos.cc (revision 4415) @@ -0,0 +1,196 @@ +/* ======================================================================== *\ +! +! * +! * 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 4415) +++ trunk/MagicSoft/Mars/mtemp/meth/MStarLocalPos.h (revision 4415) @@ -0,0 +1,86 @@ +#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