Index: /trunk/MagicSoft/Mars/mhist/MHStarMap.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHStarMap.cc (revision 3339) +++ /trunk/MagicSoft/Mars/mhist/MHStarMap.cc (revision 3340) @@ -16,7 +16,8 @@ ! ! -! Author(s): Thomas Bretz 12/2000 -! -! Copyright: MAGIC Software Development, 2000-2002 +! Author(s): Thomas Bretz 12/2000 +! Wolfgang Wittek 02/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 ! ! @@ -31,4 +32,13 @@ // from the Hillas parameters (Fill) can be enhanced. // +// For a given a shower, a series of points along its main axis are filled +// into the 2-dim histogram (x, y) of the camera plane. +// +// Before filling a point (x, y) into the histogram it is +// - shifted by (xSrc, ySrc) (the expected source position) +// - and rotated in order to compensate the rotation of the +// sky image in the camera +// +// ///////////////////////////////////////////////////////////////////////////// #include "MHStarMap.h" @@ -47,4 +57,7 @@ #include "MHillas.h" #include "MBinning.h" +#include "MMcEvt.hxx" +#include "MSrcPosCam.h" +#include "MObservatory.h" ClassImp(MHStarMap); @@ -117,4 +130,19 @@ } + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam")); + if (!fSrcPos) + { + *fLog << err << "MSrcPosCam not found... aborting" << endl; + return kFALSE; + } + + const MBinning *bins = (MBinning*)plist->FindObject("BinningStarMap"); if (!bins) @@ -134,5 +162,5 @@ if (!geom) { - *fLog << warn << dbginf << "No Camera Geometry available. Using mm-scale for histograms." << endl; + *fLog << warn << "No Camera Geometry available. Using mm-scale for histograms." << endl; return kTRUE; } @@ -141,4 +169,53 @@ return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Get the observatory location "MObservatory" from parameter list +// +Bool_t MHStarMap::ReInit(MParList *pList) +{ + MObservatory *obs = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory")); + if (!obs) + { + *fLog << err << "MObservatory not found... aborting" << endl; + return kFALSE; + } + + fSinLat = TMath::Sin(obs->GetPhi()); + fCosLat = TMath::Cos(obs->GetPhi()); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// RotationAngle +// +// calculates the angle for the rotation of the sky image in the camera; +// this angle is a function of the local coordinates (theta, phi) in rad. +// +// calculate rotation angle alpha of sky image in camera +// (see TDAS 00-11, eqs. (18) and (20)) +// +void MHStarMap::GetRotationAngle(Double_t &cosal, Double_t &sinal) +{ + const Double_t theta = fMcEvt->GetTelescopeTheta(); + const Double_t phi = fMcEvt->GetTelescopePhi(); + + const Double_t sint = TMath::Sin(theta); + const Double_t cost = TMath::Cos(theta); + + const Double_t sinp = TMath::Sin(phi); + const Double_t cosp = TMath::Cos(phi); + + const Double_t v1 = sint*sinp; + const Double_t v2 = fCosLat*cost - fSinLat*sint*cosp; + + const Double_t denom = 1./ TMath::Sqrt(v1*v1 + v2*v2); + + cosal = (fSinLat * sint) + fCosLat * cost * cosp * denom; + sinal = (fCosLat * sinp) * denom; } @@ -152,40 +229,71 @@ const MHillas &h = *(MHillas*)par; - const float delta = h.GetDelta(); - - const float m = tan(delta); - const float cosd = 1.0/sqrt(1.0+m*m); - const float sind = sqrt(1.0-cosd*cosd); - - float t = h.GetMeanY() - m*h.GetMeanX(); + const Float_t delta = h.GetDelta(); + + const Float_t m = TMath::Tan(delta); + + const Float_t cosd = 1.0/TMath::Sqrt(m*m+1); + const Float_t sind = TMath::Sqrt(1-cosd*cosd); + + Float_t t = h.GetMeanY() - m*h.GetMeanX(); + + Float_t xSource = fSrcPos->GetX(); + Float_t ySource = fSrcPos->GetY(); if (!fUseMmScale) + { t *= fMm2Deg; + xSource *= fMm2Deg; + ySource *= fMm2Deg; + } // get step size ds along the main axis of the ellipse const TAxis &axex = *fStarMap->GetXaxis(); - const float xmin = axex.GetXmin(); - const float xmax = axex.GetXmax(); + const Float_t xmin = axex.GetXmin(); + const Float_t xmax = axex.GetXmax(); // FIXME: Fixed number? - const float ds = (xmax-xmin) / 200.0; + const Float_t ds = (xmax-xmin) / 200.0; + + // Fill points along the main axis of the shower into the histogram; + // before filling + // - perform a translation by (xSource, ySource) + // - and perform a rotation to compensate the rotation of the + // sky image in the camera + Double_t cosal; + Double_t sinal; + GetRotationAngle(cosal, sinal); if (m>-1 && m<1) { - const float dx = ds * cosd; + const Float_t dx = ds * cosd; for (float x=xmin+dx/2; x<(xmax-xmin)+dx; x+=dx) - fStarMap->Fill(x, m*x+t, w); + { + const Float_t xorig = x - xSource; + const Float_t yorig = m*x+t - ySource; + + const Float_t xfill = cosal*xorig + sinal*yorig; + const Float_t yfill = -sinal*xorig + cosal*yorig; + fStarMap->Fill(xfill, yfill, w); + } } else { const TAxis &axey = *fStarMap->GetYaxis(); - const float ymin = axey.GetXmin(); - const float ymax = axey.GetXmax(); + const Float_t ymin = axey.GetXmin(); + const Float_t ymax = axey.GetXmax(); const float dy = ds * sind; for (float y=ymin+dy/2; y<(ymax-ymin)+dy; y+=dy) - fStarMap->Fill((y-t)/m, y, w); + { + const Float_t xorig = (y-t)/m - xSource; + const Float_t yorig = y - ySource; + + const Float_t xfill = cosal*xorig + sinal*yorig; + const Float_t yfill = -sinal*xorig + cosal*yorig; + fStarMap->Fill(xfill, yfill, w); + } } Index: /trunk/MagicSoft/Mars/mhist/MHStarMap.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHStarMap.h (revision 3339) +++ /trunk/MagicSoft/Mars/mhist/MHStarMap.h (revision 3340) @@ -8,9 +8,14 @@ class TH2F; class MHillas; +class MSrcPosCam; +class MMcEvt; class MHStarMap : public MH { private: - TH2F *fStarMap; //-> + MSrcPosCam *fSrcPos; //! + MMcEvt *fMcEvt; //! + + TH2F *fStarMap; //-> Float_t fMm2Deg; @@ -18,7 +23,16 @@ Bool_t fUseMmScale; + Float_t fCosLat; //! + Float_t fSinLat; //! + void PrepareDrawing() const; void Paint(Option_t *opt=""); + + void GetRotationAngle(Double_t &cosangle, Double_t &sinangle); + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + Bool_t ReInit(MParList *pList); public: @@ -29,8 +43,5 @@ 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) { return (TH1*)fStarMap; } + TH1 *GetHistByName(const TString name) { return (TH1*)fStarMap; } TH2F *GetHist() { return fStarMap; } @@ -39,5 +50,5 @@ TObject *DrawClone(Option_t *opt=NULL) const; - ClassDef(MHStarMap, 1) // Container to hold 2-dim histogram (starmap) + ClassDef(MHStarMap, 1) // Container to hold the Starmap };