Index: trunk/MagicSoft/Mars/manalysis/MHillas.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHillas.cc (revision 1426) +++ trunk/MagicSoft/Mars/manalysis/MHillas.cc (revision 1434) @@ -19,4 +19,5 @@ ! Author(s): Harald Kornmayer 1/2001 ! Author(s): Rudolf Bock 10/2001 +! Author(s): Wolfgang Wittek 6/2002 ! ! Copyright: MAGIC Software Development, 2000-2002 @@ -35,10 +36,12 @@ // Version 1: // ---------- -// fLength major axis of ellipse -// fWidth minor axis -// fDelta angle of major axis wrt x-axis -// fSize total sum of pixels -// fMeanX x of center of ellipse -// fMeanY y of center of ellipse +// fLength [mm] major axis of ellipse +// fWidth [mm] minor axis +// fDelta [rad] angle of major axis with x-axis +// by definition the major axis is pointing into +// the hemisphere x>0, thus -pi/2 < delta < pi/2 +// fSize [#CerPhot] total sum of pixels +// fMeanX [mm] x of center of ellipse +// fMeanY [mm] y of center of ellipse // // Version 2: @@ -91,9 +94,12 @@ // -------------------------------------------------------------------------- // +// Initializes the values with defaults. For the default values see the +// source code. +// void MHillas::Reset() { - fLength = 0; - fWidth = 0; - fDelta = 0; + fLength = -1; + fWidth = -1; + fDelta = 0; fSize = -1; @@ -122,9 +128,9 @@ *fLog << " - Length [mm] = " << fLength << endl; *fLog << " - Width [mm] = " << fWidth << endl; + *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl; + *fLog << " - Size [1] = " << fSize << " #CherPhot" << endl; *fLog << " - Meanx [mm] = " << fMeanX << endl; *fLog << " - Meany [mm] = " << fMeanY << endl; - *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl; *fLog << " - atg(y/x) [deg] = " << atg << endl; - *fLog << " - Size [1] = " << fSize << " #CherPhot" << endl; } @@ -207,13 +213,9 @@ // Calculate the image parameters from a Cherenkov photon event // assuming Cher.photons/pixel and pixel coordinates are given +// In case you don't call Calc from within an eventloop make sure, that +// you call the Reset member function before. // Bool_t MHillas::Calc(const MGeomCam &geom, const MCerPhotEvt &evt) { - // FIXME: MHillas::Calc is rather slow at the moment because it loops - // unnecessarily over all pixels in all its loops (s.MImgCleanStd) - // The speed could be improved very much by using Hash-Tables - // (linked lists, see THashTable and THashList, too) - // - const UInt_t npixevt = evt.GetNumPixels(); @@ -222,9 +224,17 @@ // if (npixevt <= 2) + { + *fLog << warn << "MHillas::Calc: Event has less than two pixels... skipped." << endl; return kFALSE; + } // // calculate mean value of pixel coordinates and fSize // ----------------------------------------------------- + // + // Because this are only simple sums of roughly 600 values + // with an accuracy less than three digits there should not + // be any need to use double precision (Double_t) for the + // calculation of fMeanX, fMeanY and fSize. // fMeanX = 0; @@ -235,7 +245,4 @@ fNumCorePixels = 0; - // - // FIXME! npix should be retrieved from MCerPhotEvt - // for (UInt_t i=0; i0) we + // cannot calculate Length and Width. The calculation failed + // and returnes kFALSE + // + if (corrxy==0) + { + *fLog << warn << GetDescriptor() << ": Event has CorrXY==0... skipped." << endl; + return kFALSE; + } + + // // calculate the basic Hillas parameters: orientation and size of axes // ------------------------------------------------------------------- // - const float d = (corryy - corrxx)/(corrxy*2); - - fDelta = atan(d + sqrt(d*d + 1)); - - fCosDelta = cos(fDelta); // need these in derived classes - fSinDelta = sin(fDelta); // like MHillasExt - - float axis1 = ( fCosDelta*fSinDelta*corrxy*2 + fCosDelta*fCosDelta*corrxx + fSinDelta*fSinDelta*corryy) / fSize; // [mm^2] - float axis2 = (-fCosDelta*fSinDelta*corrxy*2 + fSinDelta*fSinDelta*corrxx + fCosDelta*fCosDelta*corryy) / fSize; // [mm^2] - + // fDelta is the angle between the major axis of the ellipse and + // the x axis. It is independent of the position of the ellipse + // in the camera it has values between -pi/2 and pi/2 degrees + // + const Double_t d0 = corryy - corrxx; + const Double_t d1 = corrxy*2; + const Double_t d2 = d0 + sqrt(d0*d0 + d1*d1); + const Double_t tand = d2 / d1; + + fDelta = atan(tand); + + const Double_t s2 = tand*tand+1; + const Double_t s = sqrt(s2); + + fCosDelta = 1.0/s; // need these in derived classes + fSinDelta = tand/s; // like MHillasExt + + Double_t axis1 = (tand*tand*corryy + d2 + corrxx)/s2/fSize; + Double_t axis2 = (tand*tand*corrxx - d2 + corryy)/s2/fSize; + + // + // fLength^2 is the second moment along the major axis of the ellipse + // fWidth^2 is the second moment along the minor axis of the ellipse + // + // From the algorithm we get: fWidth <= fLength is always true + // // very small numbers can get negative by rounding - if (axis1 < 0) axis1=0; - if (axis2 < 0) axis2=0; - - fLength = sqrt(axis1); // [mm] - fWidth = sqrt(axis2); // [mm] + // + fLength = axis1<0 ? 0 : sqrt(axis1); // [mm] + fWidth = axis2<0 ? 0 : sqrt(axis2); // [mm] SetReadyToSave(); @@ -322,4 +372,5 @@ // -------------------------------------------------------------------------- // +/* void MHillas::AsciiRead(ifstream &fin) { @@ -331,5 +382,5 @@ fin >> fMeanY; } - +*/ // -------------------------------------------------------------------------- /* Index: trunk/MagicSoft/Mars/manalysis/MHillas.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHillas.h (revision 1426) +++ trunk/MagicSoft/Mars/manalysis/MHillas.h (revision 1434) @@ -15,18 +15,18 @@ private: // for description see MHillas.cc - Float_t fLength; // [mm] major axis of ellipse - Float_t fWidth; // [mm] minor axis of ellipse - Float_t fDelta; // [rad] angle of major axis with x-axis - Float_t fSize; // [#CerPhot] sum of content of all pixels (number of Cherenkov photons) - Float_t fMeanX; // [mm] x-coordinate of center of ellipse - Float_t fMeanY; // [mm] y-coordinate of center of ellipse + Float_t fLength; // [mm] major axis of ellipse + Float_t fWidth; // [mm] minor axis of ellipse + Float_t fDelta; // [rad] angle of major axis with x-axis + Float_t fSize; // [#CerPhot] sum of content of all pixels (number of Cherenkov photons) + Float_t fMeanX; // [mm] x-coordinate of center of ellipse + Float_t fMeanY; // [mm] y-coordinate of center of ellipse - Float_t fSinDelta; //! [1] sin of Delta (to be used in derived classes) - Float_t fCosDelta; //! [1] cos of Delta (to be used in derived classes) + Float_t fSinDelta; //! [1] sin of Delta (to be used in derived classes) + Float_t fCosDelta; //! [1] cos of Delta (to be used in derived classes) Short_t fNumUsedPixels; // Number of pixels which survived the image cleaning Short_t fNumCorePixels; // number of core pixels - TEllipse *fEllipse; //! Graphical Object to Display Ellipse + TEllipse *fEllipse; //! Graphical Object to Display Ellipse protected: @@ -62,5 +62,5 @@ Int_t GetNumCorePixels() const { return fNumCorePixels; } - virtual void AsciiRead(ifstream &fin); + //virtual void AsciiRead(ifstream &fin); //virtual void AsciiWrite(ofstream &fout) const; Index: trunk/MagicSoft/Mars/manalysis/MHillasSrc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHillasSrc.cc (revision 1426) +++ trunk/MagicSoft/Mars/manalysis/MHillasSrc.cc (revision 1434) @@ -19,4 +19,5 @@ ! Author(s): Harald Kornmayer 1/2001 ! Author(s): Rudolf Bock 10/2001 +! Author(s): Wolfgang Wittek 06/2002 ! ! Copyright: MAGIC Software Development, 2000-2002 @@ -32,6 +33,20 @@ // // source-dependent image parameters -// fAlpha angle between major axis and line source-to-center -// fDist distance from source to center of ellipse +// +// Version 1: +// ---------- +// fAlpha angle between major axis and line source-to-center +// fDist distance from source to center of ellipse +// +// Version 2: +// ---------- +// fHeadTail +// +// Version 3: +// ---------- +// fCosDeltaLenth cosine of angle between d and a, where +// - d is the vector from the source position to the +// center of the ellipse +// - a is a vector along the main axis of the ellipse // ///////////////////////////////////////////////////////////////////////////// @@ -57,28 +72,63 @@ } +void MHillasSrc::Reset() +{ + fDist = -1; + fAlpha = 0; + fHeadTail = 0; + fCosDeltaAlpha = 0; +} + // -------------------------------------------------------------------------- // -// calculation of source-dependent parameters +// Calculation of source-dependent parameters +// In case you don't call Calc from within an eventloop make sure, that +// you call the Reset member function before. // -void MHillasSrc::Calc(const MHillas *hillas) +Bool_t MHillasSrc::Calc(const MHillas *hillas) { - const Double_t mx = hillas->GetMeanX(); // [mm] - const Double_t my = hillas->GetMeanY(); // [mm] + fHillas = hillas; - const Double_t sx = mx - fSrcPos->GetX(); // [mm] - const Double_t sy = my - fSrcPos->GetY(); // [mm] + const Double_t mx = GetMeanX(); // [mm] + const Double_t my = GetMeanY(); // [mm] - const Double_t sd = sin(hillas->GetDelta()); // [1] - const Double_t cd = cos(hillas->GetDelta()); // [1] + const Double_t sx = mx - fSrcPos->GetX(); // [mm] + const Double_t sy = my - fSrcPos->GetY(); // [mm] + const Double_t sd = sin(GetDelta()); // [1] + const Double_t cd = cos(GetDelta()); // [1] - fHeadTail = cd*sx + sd*sy; // [mm] - fDist = sqrt(sx*sx + sy*sy); // [mm] - fAlpha = atan((cd*sy - sd*sx)/fHeadTail); // [rad] - fAlpha *= kRad2Deg; // [deg] + const Double_t tand = tan(GetDelta()); // [1] - fHillas = hillas; + fHeadTail = cd*sx + sd*sy; // [mm] + + // + // Distance from source position to center of ellipse. + // If the distance is 0 distance, Alpha is not specified. + // The calculation has failed and returnes kFALSE. + // + Double_t dist = sqrt(sx*sx + sy*sy); // [mm] + + if (dist==0) + { + *fLog << warn << GetDescriptor() << ": Event has Dist==0... skipped." << endl; + return kFALSE; + } + + fDist = dist; + // [mm] + // + // Calculate Alpha and Cosda = cos(d,a) + // The sign of Cosda will be used for quantities containing + // a head-tail information + // + const Double_t arg = (sy-tand*sx) / (fDist*sqrt(tand*tand+1)); + + fAlpha = asin(arg)*kRad2Deg; // [deg] + fCosDeltaAlpha = fHeadTail/fDist; // [1] SetReadyToSave(); + + return kTRUE; } @@ -87,7 +137,8 @@ *fLog << all; *fLog << "Source dependant Image Parameters (" << GetName() << ")" << endl; - *fLog << " - Dist = " << fDist << " mm" << endl; - *fLog << " - Alpha = " << fAlpha << " deg" << endl; - *fLog << " - HeadTail = " << fHeadTail << " mm" << endl; + *fLog << " - Dist [mm] = " << fDist << endl; + *fLog << " - Alpha [deg] = " << fAlpha << endl; + *fLog << " - HeadTail [mm] = " << fHeadTail << endl; + *fLog << " - CosDeltaAlpha = " << fCosDeltaAlpha << endl; } @@ -96,4 +147,5 @@ // overloaded MParContainer to read MHillasSrc from an ascii file // +/* void MHillasSrc::AsciiRead(ifstream &fin) { @@ -102,5 +154,5 @@ fin >> fHeadTail; } - +*/ // ----------------------------------------------------------------------- // Index: trunk/MagicSoft/Mars/manalysis/MHillasSrc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHillasSrc.h (revision 1426) +++ trunk/MagicSoft/Mars/manalysis/MHillasSrc.h (revision 1434) @@ -14,7 +14,8 @@ const MSrcPosCam *fSrcPos; //! Source position in the camere - Float_t fAlpha; // [deg] angle of major axis with vector to src - Float_t fDist; // [mm] distance from src to center of ellipse - Float_t fHeadTail; // [mm] + Float_t fAlpha; // [deg] angle of major axis with vector to src + Float_t fDist; // [mm] distance between src and center of ellipse + Float_t fHeadTail; // [mm] + Float_t fCosDeltaAlpha; // [1] cosine of angle between d and a public: @@ -24,28 +25,26 @@ const MSrcPosCam *GetSrcPos() const { return fSrcPos; } - Float_t GetLength() const { return fHillas->GetLength(); } - Float_t GetWidth() const { return fHillas->GetWidth(); } - Float_t GetDelta() const { return fHillas->GetDelta(); } - Float_t GetSize() const { return fHillas->GetSize(); } - Float_t GetMeanX() const { return fHillas->GetMeanX(); } - Float_t GetMeanY() const { return fHillas->GetMeanY(); } - Float_t GetAlpha() const { return fAlpha; } - Float_t GetDist() const { return fDist; } - Float_t GetHeadTail() const { return fHeadTail; } + void Reset(); + + Float_t GetLength() const { return fHillas->GetLength(); } + Float_t GetWidth() const { return fHillas->GetWidth(); } + Float_t GetDelta() const { return fHillas->GetDelta(); } + Float_t GetSize() const { return fHillas->GetSize(); } + Float_t GetMeanX() const { return fHillas->GetMeanX(); } + Float_t GetMeanY() const { return fHillas->GetMeanY(); } + Float_t GetAlpha() const { return fAlpha; } + Float_t GetDist() const { return fDist; } + Float_t GetHeadTail() const { return fHeadTail; } + Float_t GetCosDeltaAlpha() const { return fCosDeltaAlpha; } void Print(Option_t *opt=NULL) const; - virtual void Calc(const MHillas *hillas); + virtual Bool_t Calc(const MHillas *hillas); - virtual void AsciiRead(ifstream &fin); + //virtual void AsciiRead(ifstream &fin); //virtual void AsciiWrite(ofstream &fout) const; - ClassDef(MHillasSrc, 2) // Container to hold source position dependant parameters + ClassDef(MHillasSrc, 3) // Container to hold source position dependant parameters }; -/* - Version 1: fAlpha, fDist - Version 2: fAlpha, fDist, fSign - */ - #endif Index: trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc (revision 1426) +++ trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc (revision 1434) @@ -91,7 +91,5 @@ Bool_t MHillasSrcCalc::Process() { - fHillasSrc->Calc(fHillas); - - return kTRUE; + return fHillasSrc->Calc(fHillas) ? kTRUE : kCONTINUE; }