Index: trunk/MagicSoft/Mars/mtemp/MHTelAxisFromStars.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MHTelAxisFromStars.cc (revision 4545) +++ trunk/MagicSoft/Mars/mtemp/MHTelAxisFromStars.cc (revision 4545) @@ -0,0 +1,571 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 07/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHTelAxisFromStars +// +// This class contains histograms for MTelAxisFromStars +// +///////////////////////////////////////////////////////////////////////////// +#include "MHTelAxisFromStars.h" + +#include + +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MBinning.h" + +#include "MStarLocalCam.h" +#include "MStarLocalPos.h" +#include "MSkyCamTrans.h" +#include "MSrcPosCam.h" + + +ClassImp(MHTelAxisFromStars); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Setup the histograms +// +MHTelAxisFromStars::MHTelAxisFromStars(const char *name, const char *title) + : fMm2Deg(1), fUseMmScale(kFALSE) +{ + fName = name ? name : "MHTelAxisFromStars"; + fTitle = title ? title : "Plots for MTelAxisFromStars"; + + fNStars = new TH1D("NStars", "No. of stars", 10, -0.5, 9.5); + fNStars->SetDirectory(NULL); + fNStars->SetXTitle("Numder of stars"); + fNStars->SetYTitle("Counts"); + + fNdoF = new TH1D("NdoF", "No. of degrees of freedom", 20, -0.5, 19.5); + fNdoF->SetDirectory(NULL); + fNdoF->SetXTitle("Numder of degrees of freedom"); + fNdoF->SetYTitle("Counts"); + + fChi2Prob = new TH1D("Chi2-Prob", "Chi2 probability", 40, 0.0, 1.0); + fChi2Prob->SetDirectory(NULL); + fChi2Prob->SetXTitle("Chi2 probability"); + fChi2Prob->SetYTitle("Counts"); + + + fNumIter = new TH1D("NumIter", "Number of iterations", 50, -0.5, 49.5); + fNumIter->SetDirectory(NULL); + fNumIter->SetXTitle("Number of iterations"); + fNumIter->SetYTitle("Counts"); + + + fLambda = new TH1D("Lambda", "Scale factor lambda", 80, 0.90, 1.10); + fLambda->SetDirectory(NULL); + fLambda->SetXTitle("Scale factor lambda"); + fLambda->SetYTitle("Counts"); + + + fAlfa = new TH1D("Alfa", "Rotation angle alfa", 100, -2.5, 2.5); + fAlfa->SetDirectory(NULL); + fAlfa->SetXTitle("Rotation angle alfa [\\circ]"); + fAlfa->SetYTitle("Counts"); + + + fShift = new TH2D("Shift", "Sky-Cam transformnation : (x,y) shift", + 51, -445, 445, 51, -445, 445); + fShift->SetDirectory(NULL); + fShift->SetZTitle("Counts"); + fShift->SetXTitle("x-shift [\\circ]"); + fShift->SetYTitle("y-shift [\\circ]"); + + + fEstPos1 = new TH2D("EstPos1", "Estimated position1", + 51, -445, 445, 51, -445, 445); + fEstPos1->SetDirectory(NULL); + fEstPos1->SetZTitle("Counts"); + fEstPos1->SetXTitle("Estimated position1 x [\\circ]"); + fEstPos1->SetYTitle("Estimated position1 y [\\circ]"); + + fEstPos2 = new TH2D("EstPos2", "Estimated position2", + 51, -445, 445, 51, -445, 445); + fEstPos2->SetDirectory(NULL); + fEstPos2->SetZTitle("Counts"); + fEstPos2->SetXTitle("Estimated position2 x [\\circ]"); + fEstPos2->SetYTitle("Estimated position2 y [\\circ]"); + + // default input type : results from correlated Gauss fit + fInputType = 2; +} + +// -------------------------------------------------------------------------- +// +// Set the type of the input +// +// type = 0 calculated star positions (by averaging) +// type = 1 fitted star positions (by uncorrelated Gauss fit) +// type = 2 fitted star positions (by correlated Gauss fit) +// +void MHTelAxisFromStars::SetInputType(Int_t type) +{ + *fLog << all << "MHTelAxisFromStars::SetInputType; input type is set equal to : " + << type << endl; + + fInputType = type; +} + + +// -------------------------------------------------------------------------- +// +// Delete the histograms +// +MHTelAxisFromStars::~MHTelAxisFromStars() +{ + delete fNStars; + delete fNdoF; + delete fChi2Prob; + delete fNumIter; + delete fLambda; + delete fAlfa; + + delete fShift; + delete fEstPos1; + delete fEstPos2; +} + +// -------------------------------------------------------------------------- +// +// Setup the Binning for the histograms. +// +// Find the pointers to the containers. +// +// Use this function if you want to set the conversion factor which +// is used to convert the mm-scale in the camera plane into the deg-scale. +// The conversion factor is part of the camera geometry. Please create a +// corresponding MGeomCam container. +// +Bool_t MHTelAxisFromStars::SetupFill(const MParList *plist) +{ + + fStarLocalCam = (MStarLocalCam*)plist->FindObject("MStarLocalCam", "MStarLocalCam"); + if (!fStarLocalCam) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MStarLocalCam' not found... aborting." << endl; + return kFALSE; + } + + + fSourceLocalCam = (MStarLocalCam*)plist->FindObject("MSourceLocalCam", "MStarLocalCam"); + if (!fSourceLocalCam) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MSourceLocalCam' not found... aborting." << endl; + return kFALSE; + } + + + + fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam")); + if (!fSrcPos) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; MSrcPosCam not found... aborting" << endl; + return kFALSE; + } + + fSkyCamTrans = (MSkyCamTrans*)plist->FindObject(AddSerialNumber("MSkyCamTrans")); + if (!fSkyCamTrans) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; MSkyCamTrans not found... aborting" << endl; + return kFALSE; + } + + + //------------------------------------------------------------------ + const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam"); + if (!geom) + { + *fLog << warn << GetDescriptor() + << ": No Camera Geometry available. Using mm-scale for histograms." + << endl; + SetMmScale(kTRUE); + } + else + { + fMm2Deg = geom->GetConvMm2Deg(); + SetMmScale(kFALSE); + } + + ApplyBinning(*plist, "NStars", fNStars); + ApplyBinning(*plist, "NdoF", fNdoF); + ApplyBinning(*plist, "Chi2Prob", fChi2Prob); + ApplyBinning(*plist, "NumIter", fNumIter); + ApplyBinning(*plist, "Lambda", fLambda); + ApplyBinning(*plist, "Alfa", fAlfa); + + const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera"); + if (!bins) + { + float r = geom ? geom->GetMaxRadius() : 600; + r *= 0.9; + if (!fUseMmScale) + r *= fMm2Deg; + + MBinning b; + b.SetEdges(61, -r, r); + SetBinning(fShift, &b, &b); + SetBinning(fEstPos1, &b, &b); + SetBinning(fEstPos2, &b, &b); + } + else + { + SetBinning(fShift, bins, bins); + SetBinning(fEstPos1, bins, bins); + SetBinning(fEstPos2, bins, bins); + } + + //----------------------------------------------- + // copy names from MStarLocalCam to the histograms + MStarLocalPos *starSource = 0; + TIter nextSource(fSourceLocalCam->GetList()); + + starSource = (MStarLocalPos*)nextSource(); + if (starSource) + { + fEstPos1->SetName(starSource->GetName()); + + starSource = (MStarLocalPos*)nextSource(); + if (starSource) + fEstPos2->SetName(starSource->GetName()); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Use this function to setup your own conversion factor between degrees +// and millimeters. The conversion factor should be the one calculated in +// MGeomCam. Use this function with Caution: You could create wrong values +// by setting up your own scale factor. +// +void MHTelAxisFromStars::SetMm2Deg(Float_t mmdeg) +{ + if (mmdeg<0) + { + *fLog << warn << dbginf + << "Warning - Conversion factor < 0 - nonsense. Ignored." + << endl; + return; + } + + if (fMm2Deg>=0) + *fLog << warn << dbginf + << "Warning - Conversion factor already set. Overwriting" + << endl; + + fMm2Deg = mmdeg; +} + +// -------------------------------------------------------------------------- +// +// With this function you can convert the histogram ('on the fly') between +// degrees and millimeters. +// +void MHTelAxisFromStars::SetMmScale(Bool_t mmscale) +{ + if (fUseMmScale == mmscale) + return; + + if (fMm2Deg<0) + { + *fLog << warn << dbginf + << "Warning - Sorry, no conversion factor for conversion available." + << endl; + return; + } + + const Double_t scale = mmscale ? 1./fMm2Deg : fMm2Deg; + MH::ScaleAxis(fShift, scale, scale); + MH::ScaleAxis(fEstPos1, scale, scale); + MH::ScaleAxis(fEstPos2, scale, scale); + + if (mmscale) + { + fShift->SetXTitle("x-shift [mm]"); + fShift->SetYTitle("y-shift [mm]"); + + fEstPos1->SetXTitle("Estimated position1 x [mm]"); + fEstPos1->SetYTitle("Estimated position1 y [mm]"); + + fEstPos2->SetXTitle("Estimated position2 x [mm]"); + fEstPos2->SetYTitle("Estimated position2 y [mm]"); + } + else + { + fShift->SetXTitle("x-shift [\\circ]"); + fShift->SetYTitle("y-shift [\\circ]"); + + fEstPos1->SetXTitle("Estimated position1 x [\\circ]"); + fEstPos1->SetYTitle("Estimated position1 y [\\circ]"); + + fEstPos2->SetXTitle("Estimated position2 x [\\circ]"); + fEstPos2->SetYTitle("Estimated position2 y [\\circ]"); + } + + fUseMmScale = mmscale; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms +// +// +Bool_t MHTelAxisFromStars::Fill(const MParContainer *par, const Stat_t w) +{ + Int_t Ndof = fSkyCamTrans->GetNdof(); + if (Ndof < 0) + return kTRUE; + + const Double_t scale = fUseMmScale ? 1 : fMm2Deg; + + + fNStars ->Fill(fSkyCamTrans->GetNStars(), w); + fNdoF ->Fill(fSkyCamTrans->GetNdof(), w); + fChi2Prob ->Fill(fSkyCamTrans->GetChiSquareProb(), w); + fNumIter ->Fill(fSkyCamTrans->GetNumIter(), w); + fLambda ->Fill(fSkyCamTrans->GetLambda(), w); + fAlfa ->Fill(fSkyCamTrans->GetAlfa(), w); + + + Double_t lowx; + Double_t higx; + Double_t lowy; + Double_t higy; + Double_t x; + Double_t y; + Int_t nbix; + Int_t nbiy; + + nbix = fShift->GetXaxis()->GetNbins(); + lowx = fShift->GetXaxis()->GetBinLowEdge(1); + higx = fShift->GetXaxis()->GetBinLowEdge(nbix+1); + + nbiy = fShift->GetYaxis()->GetNbins(); + lowy = fShift->GetYaxis()->GetBinLowEdge(1); + higy = fShift->GetYaxis()->GetBinLowEdge(nbiy+1); + + x = scale * (fSkyCamTrans->GetShiftD())[0]; + y = scale * (fSkyCamTrans->GetShiftD())[1]; + if (x>lowx && xlowy && yFill(x, y, w); + } + + + MStarLocalPos *starSource = 0; + TIter nextSource(fSourceLocalCam->GetList()); + + starSource = (MStarLocalPos*)nextSource(); + if (starSource) + { + if (fInputType == 2) + { + x = scale * starSource->GetMeanXCGFit(); + y = scale * starSource->GetMeanYCGFit(); + } + else if (fInputType == 1) + { + x = scale * starSource->GetMeanXFit(); + y = scale * starSource->GetMeanYFit(); + } + else + { + x = scale * starSource->GetMeanXCalc(); + y = scale * starSource->GetMeanYCalc(); + } + + if (x>lowx && xlowy && yFill(x, y, w); + } + + + starSource = (MStarLocalPos*)nextSource(); + if (starSource) + { + if (fInputType == 2) + { + x = scale * starSource->GetMeanXCGFit(); + y = scale * starSource->GetMeanYCGFit(); + } + else if (fInputType == 1) + { + x = scale * starSource->GetMeanXFit(); + y = scale * starSource->GetMeanYFit(); + } + else + { + x = scale * starSource->GetMeanXCalc(); + y = scale * starSource->GetMeanYCalc(); + } + + if (x>lowx && xlowy && yFill(x, y, w); + } + } + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Setup an inversed deep blue sea palette. +// +void MHTelAxisFromStars::SetColors() const +{ + gStyle->SetPalette(51, NULL); + Int_t c[50]; + for (int i=0; i<50; i++) + c[49-i] = gStyle->GetColorPalette(i); + gStyle->SetPalette(50, c); +} + +// -------------------------------------------------------------------------- +// +// Creates a new canvas and draws the histograms. +// Be careful: The histograms belongs to this object and won't get deleted +// together with the canvas. +// + +void MHTelAxisFromStars::Draw(Option_t *opt) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(""); + + //TCanvas *pad = new TCanvas("TelAxisFromStars", "TelAxis plots", 900, 900); + //gROOT->SetSelectedPad(NULL); + + pad->Divide(3,3); + + pad->cd(1); + gPad->SetBorderMode(0); + fNStars->Draw(opt); + + pad->cd(2); + gPad->SetBorderMode(0); + fNdoF->Draw(opt); + + pad->cd(3); + gPad->SetBorderMode(0); + fChi2Prob->Draw(opt); + + pad->cd(4); + gPad->SetBorderMode(0); + fNumIter->Draw(opt); + + pad->cd(5); + gPad->SetBorderMode(0); + fLambda->Draw(opt); + + pad->cd(6); + gPad->SetBorderMode(0); + fAlfa->Draw(opt); + + pad->cd(7); + gPad->SetBorderMode(0); + //SetColors(); + //fShift->Draw("colz"); + fShift->Draw(""); + + pad->cd(8); + gPad->SetBorderMode(0); + //SetColors(); + fEstPos1->Draw(""); + + pad->cd(9); + gPad->SetBorderMode(0); + //SetColors(); + fEstPos2->Draw(""); + + pad->Modified(); + pad->Update(); +} + +//--------------------------------------------------------------------- +// + +TH1 *MHTelAxisFromStars::GetHistByName(const TString name) +{ + if (name.Contains("NStars", TString::kIgnoreCase)) + return fNStars; + if (name.Contains("NdoF", TString::kIgnoreCase)) + return fNdoF; + if (name.Contains("Chi2Prob", TString::kIgnoreCase)) + return fChi2Prob; + if (name.Contains("NumIter", TString::kIgnoreCase)) + return fNumIter; + if (name.Contains("Lambda", TString::kIgnoreCase)) + return fLambda; + if (name.Contains("Alfa", TString::kIgnoreCase)) + return fAlfa; + + if (name.Contains("Shift", TString::kIgnoreCase)) + return fShift; + if (name.Contains("EstPos1", TString::kIgnoreCase)) + return fEstPos1; + if (name.Contains("EstPos2", TString::kIgnoreCase)) + return fEstPos2; + + return NULL; +} + +//--------------------------------------------------------------------- +// + +void MHTelAxisFromStars::Paint(Option_t *opt) +{ + SetColors(); + MH::Paint(); +} + +//========================================================================== + + + + + + + + Index: trunk/MagicSoft/Mars/mtemp/MHTelAxisFromStars.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MHTelAxisFromStars.h (revision 4545) +++ trunk/MagicSoft/Mars/mtemp/MHTelAxisFromStars.h (revision 4545) @@ -0,0 +1,79 @@ +#ifndef MARS_MHTelAxisFromStars +#define MARS_MHTelAxisFromStars + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; +class MTelAxisFromStars; +class MSrcPosCam; +class MSkyCamTrans; +class MStarLocalCam; + +class MHTelAxisFromStars : public MH +{ +private: + + TH1D *fNStars; //-> Number of stars + TH1D *fNdoF; //-> Number of degrees of freedom + TH1D *fChi2Prob; //-> Chi2 probability + TH1D *fNumIter; //-> Number of iterations + TH1D *fLambda; //-> Scale factor lambda + TH1D *fAlfa; //-> Rotation angle alfa + + TH2D *fShift; //-> Shift between Sky and Camera system + TH2D *fEstPos1; //-> Estimated position 1 + TH2D *fEstPos2; //-> Estimated position 2 + + MStarLocalCam *fStarLocalCam; + MStarLocalCam *fSourceLocalCam; + MSrcPosCam *fSrcPos; + MSkyCamTrans *fSkyCamTrans; + + void SetColors() const; + + Float_t fMm2Deg; + Bool_t fUseMmScale; + Int_t fInputType; + + void Paint(Option_t *opt=""); + +public: + MHTelAxisFromStars(const char *name=NULL, const char *title=NULL); + ~MHTelAxisFromStars(); + + void SetInputType(Int_t type=2); + void SetMmScale(Bool_t mmscale=kTRUE); + virtual void SetMm2Deg(Float_t mmdeg); + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + + TH1 *GetHistByName(const TString name); + + TH1D *GetHistNStars() { return fNStars; } + TH1D *GetHistNdoF() { return fNdoF; } + TH1D *GetHistChi2Prob() { return fChi2Prob; } + + TH1D *GetHistNumIter() { return fNumIter; } + TH1D *GetHistLambda() { return fLambda; } + + TH1D *GetHistAlfa() { return fAlfa; } + + TH2D *GetHistShift() { return fShift; } + TH2D *GetHistEstPos1() { return fEstPos1; } + TH2D *GetHistEstPos2() { return fEstPos2; } + + void Draw(Option_t *opt=NULL); + + ClassDef(MHTelAxisFromStars, 1) // Container which holds histograms for MTelAxisFromStars +}; + +#endif + + + + + Index: trunk/MagicSoft/Mars/mtemp/MSkyCamTrans.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MSkyCamTrans.cc (revision 4545) +++ trunk/MagicSoft/Mars/mtemp/MSkyCamTrans.cc (revision 4545) @@ -0,0 +1,118 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek , 7/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////////// +// +// MSkyCamTrans +// +// container holding the parameters of the transformation from +// sky directions 'a' (projected onto the camera) +// to positions 'b' in the camera +// +// ( cos(fAlfa) -sin(fAlfa) ) +// b = fLambda * fA * a + fD fA = ( ) +// ^ ^ ^ ( sin(fAlfa) cos(fAlfa) ) +// | | | +// scale rotation shift +// factor matrix +// +// fNumIter number of iterations +// fNdof number of degrees of freedom +// fChiSquare chi-square value +// fChiSquareProb chi-square probability +// +// The units are assumed to be +// [degrees] for fAlfa +// [mm] for a, b, fD +// [1] for fLambda +// +// +// The container is filled by the task 'MTelAxisFromStars' +// +/////////////////////////////////////////////////////////////////////////// + + +#include "MSkyCamTrans.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSkyCamTrans); + +using namespace std; + +// --------------------------------------------------------------------- +// +// +// +MSkyCamTrans::MSkyCamTrans(const char *name, const char *title) +{ + fName = name ? name : "MSkyCamTrans"; + fTitle = title ? title : "Sky-Cam transformation parameters"; +} + +// --------------------------------------------------------------------- +// +// +// +void MSkyCamTrans::SetParameters(Double_t &lambda, Double_t &alfa, + Double_t a[2][2], Double_t d[2], Double_t errd[2][2], + Int_t &nstars, Int_t &numiter, + Int_t &ndof, Double_t &chisquare, Double_t &chisquareprob) +{ + fLambda = lambda; + fAlfa = alfa; + + fA[0][0] = a[0][0]; + fA[0][1] = a[0][1]; + fA[1][0] = a[1][0]; + fA[1][1] = a[1][1]; + + fD[0] = d[0]; + fD[1] = d[1]; + + fErrD[0][0] = errd[0][0]; + fErrD[0][1] = errd[0][1]; + fErrD[1][0] = errd[1][0]; + fErrD[1][1] = errd[1][1]; + + fNStars = nstars; + fNumIter = numiter; + fNdof = ndof; + fChiSquare = chisquare; + fChiSquareProb = chisquareprob; +} +// --------------------------------------------------------------------- + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mtemp/MSkyCamTrans.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MSkyCamTrans.h (revision 4545) +++ trunk/MagicSoft/Mars/mtemp/MSkyCamTrans.h (revision 4545) @@ -0,0 +1,55 @@ +#ifndef MARS_MSkyCamTrans +#define MARS_MSkyCamTrans + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MSkyCamTrans : public MParContainer +{ +private: + + //Parameters of transformation from sky to camera + + Double_t fLambda; // scale factor lambda + Double_t fAlfa; // rotation angle alfa [degrees] + Double_t fA[2][2]; // rotation matrix A + Double_t fD[2]; // shift D [mm] + Double_t fErrD[2][2]; // error matrix of shift D [mm*mm] + + Int_t fNStars; + Int_t fNumIter; + Int_t fNdof; + Double_t fChiSquare; + Double_t fChiSquareProb; + +public: + + MSkyCamTrans(const char *name=NULL, const char *title=NULL); + + void SetParameters(Double_t &, Double_t &, + Double_t[2][2], Double_t[2], Double_t[2][2], + Int_t &, Int_t &, Int_t &, Double_t &, Double_t &); + + Int_t GetNStars() { return fNStars; } + Int_t GetNumIter() { return fNumIter; } + Int_t GetNdof() { return fNdof; } + Double_t GetChiSquare() { return fChiSquare; } + Double_t GetChiSquareProb() { return fChiSquareProb; } + Double_t GetLambda() { return fLambda; } + Double_t GetAlfa() { return fAlfa; } + + Double_t *GetRotationMatrix() { return &fA[0][0]; } + Double_t *GetShiftD() { return &fD[0]; } + Double_t *GetErrMatrixShiftD() { return &fErrD[0][0]; } + + ClassDef(MSkyCamTrans, 1) // Container holding the sky-camera transformation parameters +}; + +#endif + + + + + + Index: trunk/MagicSoft/Mars/mtemp/MStarLocalPos.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MStarLocalPos.cc (revision 4533) +++ trunk/MagicSoft/Mars/mtemp/MStarLocalPos.cc (revision 4545) @@ -16,6 +16,7 @@ ! ! -! Author(s): Javier López , 4/2004 -! Author(s): Robert Wagner, 7/2004 +! Author(s): Javier López , 4/2004 +! Robert Wagner, 7/2004 +! Wolfgang Wittek, 8/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 @@ -53,5 +54,4 @@ //Info from calculation - fMagCalc = 0.; fMaxCalc = 0.; @@ -61,6 +61,5 @@ fSigmaMajorAxisCalc = 0.; - //Info from fit - + //Info from uncorrelated Gauss fit fMagFit = 0.; fMaxFit = 0.; @@ -69,6 +68,22 @@ fSigmaMinorAxisFit = 0.; fSigmaMajorAxisFit = 0.; + fChiSquare = 0.; fNdof = 1; + + //Info from correlated Gauss fit + fMagCGFit = 0.; + fMaxCGFit = 0.; + fMeanXCGFit = 0.; + fMeanYCGFit = 0.; + fSigmaXCGFit = 0.; + fSigmaYCGFit = 0.; + fCorrXYCGFit = 0.; + fXXErrCGFit = 0.; + fXYErrCGFit = 0.; + fYYErrCGFit = 0.; + + fChiSquareCGFit = 0.; + fNdofCGFit = 1; } @@ -81,5 +96,6 @@ } -void MStarLocalPos::SetCalcValues(Float_t mag, Float_t max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis) +void MStarLocalPos::SetCalcValues(Float_t mag, Float_t max, + Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis) { fMagCalc = mag; @@ -91,5 +107,7 @@ } -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) +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; @@ -103,5 +121,8 @@ } -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, Float_t xx, Float_t xy, Float_t yy) +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, + Float_t xx, Float_t xy, Float_t yy) { SetFitValues(mag, max, x, y, sigmaMinorAxis, sigmaMajorAxis, chiSquare, ndof); @@ -109,4 +130,25 @@ fYYErr = yy; fXYErr = xy; +} + +void MStarLocalPos::SetCGFitValues( + Float_t mag, Float_t max, Float_t x, Float_t y, + Float_t sigmaX, Float_t sigmaY, Float_t correlation, + Float_t xx, Float_t xy, Float_t yy, + Float_t chiSquare, Int_t ndof) +{ + fMagCGFit = mag; + fMaxCGFit = max; + fMeanXCGFit = x; + fMeanYCGFit = y; + fSigmaXCGFit = sigmaX; + fSigmaYCGFit = sigmaY; + fCorrXYCGFit = correlation; + fXXErrCGFit = xx; + fXYErrCGFit = xy; + fYYErrCGFit = yy; + + fChiSquareCGFit = chiSquare; + fNdofCGFit = ndof; } @@ -119,5 +161,4 @@ { //Print a cross in the expected position - TMarker mexp(fXExp, fYExp, 29); mexp.SetMarkerSize(3); @@ -129,9 +170,9 @@ TEllipse ecalc(fMeanXCalc, fMeanYCalc, fSigmaMinorAxisCalc, fSigmaMajorAxisCalc, 0, 360, 0); ecalc.SetLineWidth(3); - ecalc.SetLineColor(kRed); + ecalc.SetLineColor(kBlue); ecalc.Paint(); } - if (fSigmaMinorAxisFit>0. || fSigmaMajorAxisFit>0.) + if (fSigmaMinorAxisFit>0. && fSigmaMajorAxisFit>0.) { TEllipse efit(fMeanXFit, fMeanYFit, fSigmaMinorAxisFit, fSigmaMajorAxisFit, 0, 360, 0); @@ -140,4 +181,38 @@ efit.Paint(); } + + if (fSigmaXCGFit>0. && fSigmaYCGFit>0.) + { + //Print a cross in the fitted position + //TMarker mCGFit(fMeanXCGFit, fMeanYCGFit, 3); + //mCGFit.SetMarkerSize(3); + //mCGFit.SetMarkerColor(1); + //mCGFit.Paint(); + + Double_t cxx = fSigmaXCGFit*fSigmaXCGFit; + Double_t cyy = fSigmaYCGFit*fSigmaYCGFit; + Double_t d = cyy - cxx; + Double_t cxy = fCorrXYCGFit * fSigmaXCGFit * fSigmaYCGFit; + Double_t tandel; + if (cxy != 0.0) + tandel = ( d + sqrt(d*d + 4.0*cxy*cxy) ) / (2.0*cxy); + else + tandel = 0.0; + + Double_t sindel = tandel / sqrt(1.0 + tandel*tandel); + Double_t delta = TMath::ASin(sindel); + + Double_t major = (cxx + 2.0*tandel*cxy + tandel*tandel*cyy) + / (1.0 + tandel*tandel); + + Double_t minor = (tandel*tandel*cxx - 2.0*tandel*cxy + cyy) + / (1.0 + tandel*tandel); + + TEllipse efit(fMeanXCGFit, fMeanYCGFit, sqrt(major), sqrt(minor), + 0, 360, delta*kRad2Deg); + efit.SetLineWidth(3); + efit.SetLineColor(kMagenta); + efit.Paint(); + } } @@ -159,4 +234,5 @@ *fLog << inf << " Calcultated \t " << setw(4) << fMagCalc << endl; *fLog << inf << " Fitted \t " << setw(4) << fMagFit << endl; + *fLog << inf << " CGFitted \t " << setw(4) << fMagCGFit << endl; } @@ -164,6 +240,8 @@ { *fLog << inf << "Star Maximum:" << endl; - *fLog << inf << " Calcultated \t " << setw(4) << fMaxCalc << " uA" << endl; + *fLog << inf << " Calcultated \t " << setw(4) << fMaxCalc << " uA" + << endl; *fLog << inf << " Fitted \t " << setw(4) << fMaxFit << " uA" << endl; + *fLog << inf << " CGFitted \t " << setw(4) << fMaxCGFit << " uA" << endl; } @@ -171,7 +249,12 @@ { *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; + *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; + *fLog << inf << " CGFitted \t X " << setw(4) << fMeanXCGFit + << " mm \tY " << setw(4) << fMeanYCGFit << " mm" << endl; } @@ -179,6 +262,11 @@ { *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; + *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; + *fLog << inf << " CGFitted \t X " << setw(4) << fSigmaXCGFit + << " mm \tY " << setw(4) << fSigmaYCGFit << " mm \t correlation" + << setw(4) << fCorrXYCGFit << endl; } @@ -186,13 +274,30 @@ { *fLog << inf << "Star Fit Quality:" << endl; - *fLog << inf << " ChiSquare/Ndof \t " << setw(3) << fChiSquare << "/" << fNdof << endl; + *fLog << inf << " ChiSquare/Ndof \t " << setw(3) << fChiSquare + << "/" << fNdof << endl; + + *fLog << inf << "Star CGFit Quality:" << endl; + *fLog << inf << " ChiSquareCGFit/NdofCGFit \t " << setw(3) + << fChiSquareCGFit << "/" << fNdofCGFit << endl; } if (o.Contains("err", TString::kIgnoreCase) || opt == NULL) { - *fLog << inf << "Minuit Error Matrix:" << endl; - *fLog << inf << " xx,xy,yy \t " << setw(3) << fXXErr << ", " << fXYErr << ", " << fYYErr << endl; - } - - -} + *fLog << inf << "CGFit Error Matrix of (fMeanXCGFit,fMeanYCGFit) :" + << endl; + *fLog << inf << " xx,xy,yy \t " << setw(3) << fXXErrCGFit << ", " + << fXYErrCGFit << ", " << fYYErrCGFit << endl; + } + + +} +//-------------------------------------------------------------------------- + + + + + + + + + Index: trunk/MagicSoft/Mars/mtemp/MStarLocalPos.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MStarLocalPos.h (revision 4533) +++ trunk/MagicSoft/Mars/mtemp/MStarLocalPos.h (revision 4545) @@ -11,5 +11,4 @@ //Expected position on camera - Float_t fMagExp; Float_t fXExp; //[mm] @@ -17,5 +16,4 @@ //Info from calculation - Float_t fMagCalc; Float_t fMaxCalc; //[uA] @@ -25,6 +23,5 @@ Float_t fSigmaMajorAxisCalc; //[mm] - //Info from fit - + //Info from uncorrelated Gauss fit Float_t fMagFit; Float_t fMaxFit; //[uA] @@ -33,9 +30,26 @@ Float_t fSigmaMinorAxisFit; //[mm] Float_t fSigmaMajorAxisFit; //[mm] - Float_t fChiSquare; - Float_t fXXErr; //minuit error matrix elements + Float_t fXXErr; Float_t fXYErr; Float_t fYYErr; + + Float_t fChiSquare; Int_t fNdof; + + //Info from correlated Gauss fit + Float_t fMagCGFit; + Float_t fMaxCGFit; //[uA] + Float_t fMeanXCGFit; //[mm] + Float_t fMeanYCGFit; //[mm] + Float_t fSigmaXCGFit; //[mm] + Float_t fSigmaYCGFit; //[mm] + Float_t fCorrXYCGFit; // correlation coefficient + Float_t fXXErrCGFit; // error matrix of (fMeanXCGFit,fMeanYCGFit) + Float_t fXYErrCGFit; + Float_t fYYErrCGFit; + + Float_t fChiSquareCGFit; + Int_t fNdofCGFit; + public: @@ -70,14 +84,40 @@ Float_t GetSigmaMajorAxis() {return fSigmaMajorAxisFit!=0?fSigmaMajorAxisFit:fSigmaMajorAxisCalc;} - Float_t GetXXErr() {return fXXErr;} - Float_t GetXYErr() {return fXYErr;} - Float_t GetYYErr() {return fYYErr;} + // getters for the correlated Gauss fit + Float_t GetMagCGFit() {return fMagCGFit;} + Float_t GetMaxCGFit() {return fMaxCGFit;} + Float_t GetMeanXCGFit() {return fMeanXCGFit;} + Float_t GetMeanYCGFit() {return fMeanYCGFit;} + Float_t GetSigmaXCGFit() {return fSigmaXCGFit;} + Float_t GetSigmaYCGFit() {return fSigmaYCGFit;} + Float_t GetCorrXYCGFit() {return fCorrXYCGFit;} + Float_t GetXXErrCGFit() {return fXXErrCGFit;} + Float_t GetXYErrCGFit() {return fXYErrCGFit;} + Float_t GetYYErrCGFit() {return fYYErrCGFit;} + Float_t GetChiSquareCGFit() {return fChiSquareCGFit;} + UInt_t GetNdofCGFit() {return fNdofCGFit;} + Float_t GetChiSquareNdofCGFit() {return fChiSquareCGFit/fNdofCGFit;} + 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 max, Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, Float_t chi, Int_t ndof, Float_t xx, Float_t xy, Float_t yy); + + 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 max, Float_t x, Float_t y, + Float_t sigmaX, Float_t sigmaY, + Float_t chi, Int_t ndof, + Float_t xx, Float_t xy, Float_t yy); + + void SetCGFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, + Float_t sigmaX, Float_t sigmaY, Float_t correlation, + Float_t xx, Float_t xy, Float_t yy, + Float_t chi, Int_t ndof); void Paint(Option_t *opt=NULL); @@ -88,2 +128,4 @@ #endif + + Index: trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc (revision 4545) +++ trunk/MagicSoft/Mars/mtemp/MTelAxisFromStars.cc (revision 4545) @@ -0,0 +1,1069 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 07/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MTelAxisFromStars +// +// This task +// - determines the transformation from expected positions of stars +// in the camera to measured positions of these stars in the camera +// - applies this transformation to expected positions of other objects +// to obtain the estimated positions of these objects in the camera +// - puts the estimated positions into the relevant containers +// +// Input Containers : +// MStarLocalCam[MStarLocalCam], MStarLocalCamSource[MStarLocalCam] +// +// Output Containers : +// MSkyCamTrans, MSrcPosCam +// +///////////////////////////////////////////////////////////////////////////// +#include +#include + +#include + +#include "MTelAxisFromStars.h" + +#include "MParList.h" +#include "MSrcPosCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MStarLocalCam.h" +#include "MStarLocalPos.h" +#include "MSkyCamTrans.h" + +ClassImp(MTelAxisFromStars); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Constructor +// +MTelAxisFromStars::MTelAxisFromStars(const char *name, const char *title) +{ + fName = name ? name : "MTelAxisFromStars"; + fTitle = title ? title : "Calculate source position from star positions"; + + // if scale factor fLambda should NOT be fixed set fFixdScaleFactor to + // -1.0; otherwise set it to the value requested + fFixedScaleFactor = 1.0; + + // if rotation angle fAlfa should NOT be fixed set fFixdRotationAngle to + // -1.0; otherwise set it to the requested value + fFixedRotationAngle = 0.0; + + // default type of input is : the result of the correlated Gauss fit + // type 0 : result from the weighted average + // type 1 : result from the uncorrelated Gauss fit + fInputType = 2; +} + +// -------------------------------------------------------------------------- +// +// Destructor +// +MTelAxisFromStars::~MTelAxisFromStars() +{ +} + +// -------------------------------------------------------------------------- +// +// Set links to containers +// + +Int_t MTelAxisFromStars::PreProcess(MParList *pList) +{ + + fStarLocalCam = (MStarLocalCam*)pList->FindObject("MStarLocalCam", "MStarLocalCam"); + if (!fStarLocalCam) + { + *fLog << err << "MTelAxisFromStars::PreProcess; container 'MStarLocalCam' not found... aborting." << endl; + return kFALSE; + } + + + fSourceLocalCam = (MStarLocalCam*)pList->FindObject("MSourceLocalCam", "MStarLocalCam"); + if (!fSourceLocalCam) + { + *fLog << err << "MTelAxisFromStars::PreProcess; container 'MSourceLocalCam' not found... continue " << endl; + } + + + fSrcPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam")); + if (!fSrcPos) + { + *fLog << err << "MTelAxisFromStars::PreProcess; MSrcPosCam not found... aborting" << endl; + return kFALSE; + } + + fSkyCamTrans = (MSkyCamTrans*)pList->FindCreateObj(AddSerialNumber("MSkyCamTrans")); + if (!fSkyCamTrans) + { + *fLog << err << "MTelAxisFromStars::PreProcess; MSkyCamTrans not found... aborting" << endl; + return kFALSE; + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set the type of the input +// +// type = 0 calculated star positions (by averaging) +// type = 1 fitted star positions (by uncorrelated Gauss fit) +// type = 2 fitted star positions (by correlated Gauss fit) +// +void MTelAxisFromStars::SetInputType(Int_t type) +{ + *fLog << all << "MTelAxisFromStars::SetInputType; input type is set equal to : " + << type << endl; + + fInputType = type; +} + +// -------------------------------------------------------------------------- +// +// Fix the scale factor fLambda +// +// +void MTelAxisFromStars::FixScaleFactorAt(Double_t lambda) +{ + *fLog << all << "MTelAxisFromStars::FixScaleFactorAt; scale factor will be fixed at : " + << lambda << endl; + + fFixedScaleFactor = lambda; +} + + +// -------------------------------------------------------------------------- +// +// Fix rotation angle fAlfa +// +// +void MTelAxisFromStars::FixRotationAngleAt(Double_t alfa) +{ + *fLog << all << "MTelAxisFromStars::FixRotationAngleAt; rotation angle will be fixed at : " + << alfa << endl; + + fFixedRotationAngle = alfa; // [degrees] +} + + +// -------------------------------------------------------------------------- +// +// Process +// +// call FindSkyCamTrans to find the Sky-Camera transformation +// call TransSkyCam to transform some sky directions +// into the camera system +// put the estimated source position into MSrcPosCam +// + +Int_t MTelAxisFromStars::Process() +{ + //Int_t run = fRun->GetRunNumber(); + //*fLog << "MTelAxisFromStars::Process; run = " << run << endl; + + //-------------------------------------- + // Define the input for FindSkyCamTrans + // + + // get the expected (axy[0], axy[1]) and the measured positions + // (bxy[0], bxy[1]) of stars in the camera from MStarLocalCam + Int_t fNumStars = fStarLocalCam->GetNumStars(); + + if (fNumStars <= 0) + return kTRUE; + + TArrayD axy[2]; + axy[0].Set(fNumStars); + axy[1].Set(fNumStars); + + TArrayD bxy[2]; + bxy[0].Set(fNumStars); + bxy[1].Set(fNumStars); + + // error matrix of bxy + TArrayD exy[2][2]; + exy[0][0].Set(fNumStars); + exy[0][1].Set(fNumStars); + exy[1][0].Set(fNumStars); + exy[1][1].Set(fNumStars); + + // transformation parameters + Double_t fLambda; + Double_t fAlfa; + Double_t fA[2][2]; + Double_t fD[2]; + Double_t fErrD[2][2]; + Int_t fNumIter; + Int_t fNdof; + Double_t fChi2; + Double_t fChi2Prob; + + MStarLocalPos *star = 0; + TIter next(fStarLocalCam->GetList()); + Int_t ix = 0; + + // loop over all stars + while ( (star = (MStarLocalPos*)next()) ) + { + axy[0][ix] = star->GetXExp(); + axy[1][ix] = star->GetYExp(); + + if (fInputType == 0) + { + // values from averaging + bxy[0][ix] = star->GetMeanXCalc(); + bxy[1][ix] = star->GetMeanYCalc(); + exy[0][0][ix]= star->GetSigmaMinorAxisCalc()*star->GetSigmaMinorAxisCalc(); + exy[0][1][ix] = 0.0; + exy[1][0][ix] = 0.0; + exy[1][1][ix]= star->GetSigmaMajorAxisCalc()*star->GetSigmaMajorAxisCalc(); + } + + else if (fInputType == 1) + { + // values from uncorrelatd Gauss fit + bxy[0][ix] = star->GetMeanXFit(); + bxy[1][ix] = star->GetMeanYFit(); + exy[0][0][ix]= star->GetSigmaMinorAxisFit()*star->GetSigmaMinorAxisFit(); + exy[0][1][ix] = 0.0; + exy[1][0][ix] = 0.0; + exy[1][1][ix]= star->GetSigmaMajorAxisFit()*star->GetSigmaMajorAxisFit(); + } + + else if (fInputType == 2) + { + // values from correlatd Gauss fit + bxy[0][ix] = star->GetMeanXCGFit(); + bxy[1][ix] = star->GetMeanYCGFit(); + + // this is the error matrix for (MeanXCGFit, MeanYCGFit); + // this is the error matrix which should be used + exy[0][0][ix] = star->GetXXErrCGFit(); + exy[0][1][ix] = star->GetXYErrCGFit(); + exy[1][0][ix] = star->GetXYErrCGFit(); + exy[1][1][ix] = star->GetYYErrCGFit(); + + // this is the error matrix constructed from SigmaXCGFit and SigmaYCGFit; + // it is used because the errors above are too small, at present + //exy[0][0][ix] = star->GetSigmaXCGFit() * star->GetSigmaXCGFit(); + //exy[0][1][ix] = star->GetCorrXYCGFit() * + // star->GetSigmaXCGFit() * star->GetSigmaYCGFit(); + //exy[1][0][ix] = exy[0][1][ix]; + //exy[1][1][ix] = star->GetSigmaYCGFit() * star->GetSigmaYCGFit(); + } + + else + { + *fLog << err << "MTelAxisFromStars::Process; type of input is not defined" + << endl; + return kFALSE; + } + + // don't include stars with undefined error + Double_t deter = exy[0][0][ix]*exy[1][1][ix] + - exy[0][1][ix]*exy[1][0][ix]; + + //*fLog << "ix ,deter, xx, xy, yy = " << ix << ": " + // << deter << ", " << exy[0][0][ix] << ", " + // << exy[0][1][ix] << ", " << exy[1][1][ix] << endl; + if (deter <= 0.0) + continue; + + //*fLog << "MTelAxisFromStars : " << endl; + //*fLog << " ix, XExp, YExp, XFit, YFit, SigmaX2, SigmaXY, SigmaY2 = " + // << ix << " : " + // << axy[0][ix] << ", " << axy[1][ix] << ", " + // << bxy[0][ix] << ", " << bxy[1][ix] << ", " + // << exy[0][0][ix] << ", " << exy[0][1][ix] << ", " + // << exy[1][1][ix] << endl; + + ix++; + } + + //-------------------------------------- + // Find the transformation from expected positions (axy[1], axy[2]) + // to measured positions (bxy[1], bxy[2]) in the camera + + Int_t fNStars = ix; + + if (ix < fNumStars) + { + // reset the sizes of the arrays + Int_t fNStars = ix; + axy[0].Set(fNStars); + axy[1].Set(fNStars); + + bxy[0].Set(fNStars); + bxy[1].Set(fNStars); + + exy[0][0].Set(fNStars); + exy[0][1].Set(fNStars); + exy[1][0].Set(fNStars); + exy[1][1].Set(fNStars); + } + + Bool_t fitOK; + if (fNStars < 1) + { + *fLog << "MTelAxisFromStars::Process; no star for MTelAxisFromStars" + << endl; + fitOK = kFALSE; + } + else + { + fitOK = FindSkyCamTrans(axy, bxy, exy, + fFixedRotationAngle, fFixedScaleFactor, fLambda, + fAlfa , fA, fD, fErrD, + fNumIter, fNdof, fChi2, fChi2Prob); + } + + if (!fitOK) + { + *fLog << err + << "MTelAxisFromStars::Process; Fit to find transformation from star to camera system failed" + << endl; + + if (fNStars > 0) + { + *fLog << err + << " fNumIter, fNdof, fChi2, fChi2Prob = " << fNumIter + << ", " << fNdof << ", " << fChi2 << ", " << fChi2Prob << endl; + } + + return kTRUE; + } + + + //-------------------------------------- + // Put the transformation parameters into the MSkyCamTrans container + + fSkyCamTrans->SetParameters(fLambda, fAlfa, fA, fD, fErrD, + fNumStars, fNumIter, fNdof, fChi2, fChi2Prob); + fSkyCamTrans->SetReadyToSave(); + + + //-------------------------------------- + // Put the estimated position, obtained by transforming the expected + // position (0, 0), into SrcPosCam + + fSrcPos->SetXY(fD[0], fD[1]); + fSrcPos->SetReadyToSave(); + + + //-------------------------------------- + // Apply the transformation to some expected positions (asxy[1], asxy[2]) + // to obtain estimated positions (bsxy[1], bsxy[2]) in the camera + // and their error matrices esxy[2][2] + + // get the expected positions (asxy[1], asxy[2]) from another MStarLocalCam + // container (with the name "MSourceLocalCam") + Int_t fNumStarsSource = fSourceLocalCam->GetNumStars(); + + //*fLog << "MTelAxisFromStars::Process; fNumStarsSource = " + // << fNumStarsSource << endl; + + if (fNumStarsSource > 0) + { + TArrayD asxy[2]; + asxy[0].Set(fNumStarsSource); + asxy[1].Set(fNumStarsSource); + + TArrayD bsxy[2]; + bsxy[0].Set(fNumStarsSource); + bsxy[1].Set(fNumStarsSource); + + TArrayD esxy[2][2]; + esxy[0][0].Set(fNumStarsSource); + esxy[0][1].Set(fNumStarsSource); + esxy[1][0].Set(fNumStarsSource); + esxy[1][1].Set(fNumStarsSource); + + MStarLocalPos *starSource = 0; + TIter nextSource(fSourceLocalCam->GetList()); + ix = 0; + while ( (starSource = (MStarLocalPos*)nextSource()) ) + { + asxy[0][ix] = starSource->GetXExp(); + asxy[1][ix] = starSource->GetYExp(); + + ix++; + } + + TransSkyCam(fLambda, fA, fD, fErrD, asxy, bsxy, esxy); + + // put the estimated positions into the MStarLocalPos container + TIter setnextSource(fSourceLocalCam->GetList()); + ix = 0; + while ( (starSource = (MStarLocalPos*)setnextSource()) ) + { + if (fInputType == 1) + { + starSource->SetFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix], + sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), fChi2, fNdof); + } + else if (fInputType == 2) + { + Double_t corr = esxy[0][1][ix]/ sqrt( esxy[0][0][ix] * esxy[1][1][ix] ); + starSource->SetCGFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix], + sqrt(esxy[0][0][ix]), sqrt(esxy[1][1][ix]), corr, + esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix], + fChi2, fNdof); + } + + ix++; + } + + } + + //-------------------------------------- + + return kTRUE; +} + +//--------------------------------------------------------------------------- +// +// FindSkyCamTrans +// +// This routine determines the transformation +// +// ( cos(alfa) -sin(alfa) ) +// b = lambda * A * a + d A = ( ) +// ^ ^ ^ ( sin(alfa) cos(alfa) ) +// | | | +// scale rotation shift +// factor matrix +// +// from sky coordinates 'a' (projected onto the camera) to camera +// coordinates 'b', using the positions of known stars in the camera. +// The latter positions may have been determined by analysing the +// DC currents in the different pixels. +// +// Input : a[2] x and y coordinates of stars projected onto the camera; +// they were obtained from (RA, dec) of the stars and +// (ThetaTel, PhiTel) and the time of observation; +// these are the 'expected positions' of stars in the camera +// b[2] 'measured positions' of these stars in the camera; +// they may have been obtained from the DC currents +// e[2][2] error matrix of b[2] +// fixedrotationangle value [in degrees] at which rotation angle +// alfa should be fixed; -1 means don't fix +// fixedscalefactor value at which scale factor lambda +// should be fixed; -1 means don't fix +// +// Output : lambda, alfadeg, A[2][2], d[2] fit results; +// parameters describing the transformation +// from 'expected positions' to the 'measured +// positions' in the camera +// errd[2][2] error matrix of d[2] +// fNumIter number of iterations +// fNdoF number of degrees of freedom +// fChi2 chi-square value +// fChi2Prob chi-square probability +// +// The units are assumed to be +// [degrees] for alfadeg +// [mm] for a, b, d +// [1] for lambda + +Bool_t MTelAxisFromStars::FindSkyCamTrans( + TArrayD a[2], TArrayD b[2], TArrayD e[2][2], + Double_t &fixedrotationang, Double_t &fixedscalefac, Double_t &lambda, + Double_t &alfadeg, Double_t A[2][2], Double_t d[2], Double_t errd[2][2], + Int_t &fNumIter, Int_t &fNdof, Double_t &fChi2, Double_t &fChi2Prob) +{ + Int_t fNumStars = a[0].GetSize(); + + //*fLog << "MTelAxisFromStars::FindSkyCamTrans; expected and measured positions :" + // << endl; + for (Int_t ix=0; ix 1.1 ) + lambda = 1.1; + + dlambda = lambda - lambdaold; + lambdaold = lambda; + + //*fLog << "num, denom, lambda, dlambda = " << num + // << ", " << denom << ", " << lambda << ", " + // << dlambda << endl; + } + + } + //------- end of iteration ***************************************** + + alfadeg = alfa * kRad2Deg; + + // calculate error matrix of d[2] + errd[0][0] = sumcinv[0][0]; + errd[0][1] = sumcinv[0][1]; + errd[1][0] = sumcinv[1][0]; + errd[1][1] = sumcinv[1][1]; + + // evaluate quality of fit + + // calculate chi2 + for (Int_t ix=0; ix +#endif + +class MSrcPosCam; +class MStarLocalCam; +class MSkyCamTrans; + +class MTelAxisFromStars : public MTask +{ + private: + + const MStarLocalCam *fStarLocalCam; //! + const MStarLocalCam *fSourceLocalCam; //! + + MSrcPosCam *fSrcPos; //! + MSkyCamTrans *fSkyCamTrans; //! + + Double_t fFixedScaleFactor; //! + Double_t fFixedRotationAngle; //! [degrees] + + Int_t fInputType; //! type of input + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + public: + + MTelAxisFromStars(const char *name=NULL, const char *title=NULL); + ~MTelAxisFromStars(); + + void FixScaleFactorAt(Double_t lambda = 1.0); + void FixRotationAngleAt(Double_t alfa = 0.0); // alfa in [degrees] + + void SetInputType(Int_t type = 2); + + Bool_t FindSkyCamTrans(TArrayD[2], TArrayD[2], TArrayD[2][2], + Double_t &, Double_t &, Double_t &, + Double_t &, Double_t[2][2] , Double_t[2], Double_t[2][2], + Int_t &, Int_t &, Double_t &, Double_t &); + + void TransSkyCam(Double_t &, Double_t[2][2], Double_t[2], Double_t[2][2], + TArrayD[2], TArrayD[2], TArrayD[2][2]); + + ClassDef(MTelAxisFromStars, 0) // Task to calculate the source position from star positions +}; + +#endif + + + + + + + + Index: trunk/MagicSoft/Mars/mtemp/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/Makefile (revision 4533) +++ trunk/MagicSoft/Mars/mtemp/Makefile (revision 4545) @@ -40,5 +40,9 @@ MStarLocalPos.cc \ MStarLocalCam.cc \ - MFindStars.cc + MFindStars.cc \ + MTelAxisFromStars.cc \ + MHTelAxisFromStars.cc \ + MSkyCamTrans.cc \ + MSourceDirections.cc ############################################################ Index: trunk/MagicSoft/Mars/mtemp/TempLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/TempLinkDef.h (revision 4533) +++ trunk/MagicSoft/Mars/mtemp/TempLinkDef.h (revision 4545) @@ -8,4 +8,8 @@ #pragma link C++ class MStarLocalCam+; #pragma link C++ class MFindStars+; +#pragma link C++ class MTelAxisFromStars+; +#pragma link C++ class MHTelAxisFromStars+; +#pragma link C++ class MSkyCamTrans+; +#pragma link C++ class MSourceDirections+; #endif Index: trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C (revision 4545) +++ trunk/MagicSoft/Mars/mtemp/findTelAxisFromStars.C (revision 4545) @@ -0,0 +1,265 @@ +void ReadSetup(TString fname, MAstroCamera &cam) +{ + MMcConfigRunHeader *config=0; + MGeomCam *geom=0; + + TFile file(fname); + TTree *tree = (TTree*)file.Get("RunHeaders"); + tree->SetBranchAddress("MMcConfigRunHeader", &config); + if (tree->GetBranch("MGeomCam")) + tree->SetBranchAddress("MGeomCam", &geom); + tree->GetEntry(0); + + cam.SetMirrors(*config->GetMirrors()); + cam.SetGeom(*geom); +} + +void findTelAxisFromStars(const TString filename="20040422_23213_D_Mrk421_E.root", const TString directory="/.magic/magicserv01/MAGIC/rootdata/2004_04_22/", const UInt_t numEvents = 0) +{ + gLog.SetNoColors(); + + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + MGeomCamMagic geomcam; + MCameraDC dccam; + MStarLocalCam starcam; + + //$$$$$$$$$$$$$$$$ ww + MHTelAxisFromStars htelaxis; + plist.AddToList(&htelaxis); + //$$$$$$$$$$$$$$$$ ww + + + plist.AddToList(&geomcam); + plist.AddToList(&dccam); + plist.AddToList(&starcam); + + // Reads the trees of the root file and the analysed branches + MReadReports read; + read.AddTree("Currents"); + read.AddTree("Drive"); // If you do not include Drive info, the MFindStars class + // will not use the star catalog method + read.AddFile(directory+filename); // after the reading of the trees!!! + read.AddToBranchList("MReportCurrents.*"); + read.AddToBranchList("MReportDrive.*"); + + MGeomApply geomapl; + // TString continuoslightfile = + // "/data/MAGIC/Period016/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root"; + + TString continuoslightfile = +"/.magic/magicserv01/MAGIC/rootdata/2004_04_16/20040416_22368_P_Off3c279-2CL100_E.root"; + + + Float_t mindc = 0.9; //[uA] + + MCalibrateDC dccal; + dccal.SetFileName(continuoslightfile); + dccal.SetMinDCAllowed(mindc); + + const Int_t numblind = 5; + const Short_t x[numblind] = { 47, 124, 470, 475, 571}; + const TArrayS blindpixels(numblind,(Short_t*)x); + Float_t ringinterest = 100; //[mm] + Float_t tailcut = 2.5; + UInt_t integratedevents = 1; + + // We need the MAGIC mirror geometry from a MC header: + //TString geometryfile = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root"; + TString geometryfile = +"/.magic/magicserv01/MAGIC/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root"; + + + + + // We need the Bright Star Catalog: + //TString catalogfile = "/home/rwagner/bsc5.dat"; + TString catalogfile = "mtemp/mmpi/macros/bsc5.dat"; + + //$$$$$$$$$$$$$$$$ ww + MSourceDirections sdirs; + sdirs.SetGeometryFile(geometryfile); + const Double_t ra = MAstro::Hms2Rad(11, 4, 26); + const Double_t dec = MAstro::Dms2Rad(38, 12, 36); + sdirs.SetRaDec(ra,dec); + sdirs.AddDirection(ra,dec,1,"Mkn 421"); + const Double_t ra = MAstro::Hms2Rad(11, 4, 31); + const Double_t dec = MAstro::Dms2Rad(38, 14, 29); + sdirs.AddDirection(ra,dec,1,"My_UMa 51"); + sdirs.SetRadiusFOV(3); + //$$$$$$$$$$$$$$$$ ww + + MFindStars findstars; + findstars.SetBlindPixels(blindpixels); + findstars.SetRingInterest(ringinterest); + findstars.SetDCTailCut(tailcut); + findstars.SetNumIntegratedEvents(integratedevents); + findstars.SetMinuitPrintOutLevel(-1); + findstars.SetGeometryFile(geometryfile); + findstars.SetBSCFile(catalogfile); + const Double_t ra = MAstro::Hms2Rad(11, 4, 26); + const Double_t dec = MAstro::Dms2Rad(38, 12, 36); + findstars.SetRaDec(ra,dec); + findstars.SetLimMag(8); + findstars.SetRadiusFOV(1.5); + + + + //$$$$$$$$$$$$$$$$ ww + //Int_t InputType = 0; + //findstars.SetUseCorrelatedGauss(kFALSE); + + MTelAxisFromStars telaxis; + //telaxis.FixRotationAngleAt(-1.0); + //telaxis.FixScaleFactorAt(-1.0); + //telaxis.SetInputType(InputType); + + MFillH fillhisto("MHTelAxisFromStars[MHTelAxisFromStars]",""); + //htelaxis.SetInputType(InputType); + //$$$$$$$$$$$$$$$$ ww + + + tlist.AddToList(&geomapl); + tlist.AddToList(&read); + tlist.AddToList(&dccal); + tlist.AddToList(&findstars, "Currents"); + + //$$$$$$$$$$$$$$$$ ww + tlist.AddToList(&sdirs); + tlist.AddToList(&telaxis); + tlist.AddToList(&fillhisto); + //$$$$$$$$$$$$$$$$ ww + + + + // The following lines you only need if in addition you want to display + // independent MAstroCamera output + // + // TString fname = "/mcdata/standard/camera/NSB_013/Gamma/Gamma_zbin9_90_7_1480to1489_w0.root"; + // MObservatory magic1; + // const Double_t ra = MAstro::Hms2Rad(11, 4, 26); //Mkn421 + // const Double_t dec = MAstro::Dms2Rad(38, 12, 36); + // + // MAstroCamera stars; + // ReadSetup(fname, stars); + // stars.SetLimMag(9); + // stars.SetRadiusFOV(3); + // stars.SetRaDec(ra, dec); + // stars.ReadBSC("/home/rwagner/bsc5.dat"); + // stars.SetObservatory(magic1); + + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + if (!evtloop.PreProcess()) + return; + + MHCamera display(geomcam); + display.SetPrettyPalette(); + display.Draw(); + gPad->cd(1); + starcam.Draw(); + + UInt_t numevents=0; + + while (tlist.Process()) + { + starcam.Print("namepossizchierr"); + + numevents++; + if (numevents%integratedevents==0) + { + display.SetCamContent(findstars.GetDisplay()); + gPad->Modified(); + gPad->Update(); + // This line prints the results: + // starcam.Print(); + // This is how to access the TList of stars: + // TList* starlist = starcam.GetList(); + + // This is how to iterate over stars found: + // TIter Next(starlist); + // MStarLocalPos* star; + // UInt_t starnum = 0; + // cout << filename << " "; + // cout << "Iterating over list" << endl; + // while ((star=(MStarLocalPos*)Next())) + // { + // cout << "star[" << starnum << "] "; + // cout << star->GetMeanX() << " " + // << star->GetMeanY() << " "; + // starnum++; + // } + // cout << endl; + + }//integratedevents + + MTime time; + time.Set(2004, 4, 22, 21, 51, 15); + + //superimpose star picture + // stars.SetTime(time); + // TObject *o = stars.Clone(); + // o->SetBit(kCanDelete); + // o->Draw(); + + // wait after each event + if (!HandleInput()) + break; + + } + + + + evtloop.PostProcess(); + tlist.PrintStatistics(); + + //$$$$$$$$$$$$$$$$ ww + /* + gLog << "Event loop finished; call DrawClone of MHTelAxisFromStars" << endl; + TObject *obj = plist.FindObject("MHTelAxisFromStars"); + if (obj) + { + obj->Print(); + obj->Dump(); + obj->DrawClone(); + + gLog << "Draw a clone of the histogram" << endl; + TObject *o = obj->Clone(); + o->SetBit(kCanDelete); + o->Draw(); + } + else + gLog << "address of MHTelAxisFromStars container is zero" << endl; + */ + //$$$$$$$$$$$$$$$$ ww + + + +} + + +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; +}