/* ======================================================================== *\ ! ! * ! * 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 12/2000 ! Author(s): Harald Kornmayer 1/2001 ! Author(s): Rudolf Bock 10/2001 ! Author(s): Wolfgang Wittek 6/2002 ! ! Copyright: MAGIC Software Development, 2000-2008 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MHillas // // Storage Container for image parameters // // basic image parameters // // Version 1: // ---------- // 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: // ---------- // fNumCorePixels number of pixels called core // fNumUsedPixels number of pixels which survived the cleaning // // Version 3: // ---------- // fNumCorePixels moved to MNewImagePar // fNumUsedPixels moved to MNewImagePar // fCosDelta added // fSinDelte added // ///////////////////////////////////////////////////////////////////////////// #include "MHillas.h" #include #include #include #include #include "MGeom.h" #include "MGeomCam.h" #include "MSignalPix.h" #include "MSignalCam.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHillas); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. // MHillas::MHillas(const char *name, const char *title) { fName = name ? name : "MHillas"; fTitle = title ? title : "Storage container for image parameters of one event"; Reset(); } // -------------------------------------------------------------------------- // // Initializes the values with defaults. For the default values see the // source code. // void MHillas::Reset() { fLength = -1; fWidth = -1; fDelta = 0; fSize = -1; fMeanX = 0; fMeanY = 0; } // -------------------------------------------------------------------------- // // return distance to center (coordinate origin) // Double_t MHillas::GetDist0() const { return TMath::Hypot(fMeanX, fMeanY); } // -------------------------------------------------------------------------- // // return the mean position as TVector2(meanx, meany) // TVector2 MHillas::GetMean() const { return TVector2(fMeanX, fMeanY); } // -------------------------------------------------------------------------- // // return the direction of the main axis (normalized). The main axis // always has the direction --> x>0 // TVector2 MHillas::GetNormAxis() const { return TVector2(fCosDelta, fSinDelta); } // -------------------------------------------------------------------------- // // Analytical estimation of length of border line: // U = pi*(a+b - sqrt(a*b)) // // GetBorderLine() [mm] // Double_t MHillas::GetBorderLine() const { const Double_t a = fWidth+fLength; const Double_t s = fWidth*fLength; return TMath::Pi()*(1.5*a - TMath::Sqrt(s)); } // -------------------------------------------------------------------------- // // Analytical estimation of length of border line: // A = pi*a*b // // GetArea() [mm^2] // Double_t MHillas::GetArea() const { return TMath::Pi()*fWidth*fLength; } // -------------------------------------------------------------------------- // // Print the hillas Parameters to *fLog // void MHillas::Print(Option_t *) const { const Double_t atg = atan2(fMeanX, fMeanY)*kRad2Deg; *fLog << all; *fLog << GetDescriptor() << endl; *fLog << " - Length [mm] = " << fLength << endl; *fLog << " - Width [mm] = " << fWidth << endl; *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl; *fLog << " - Size [phe] = " << fSize << endl; *fLog << " - Meanx [mm] = " << fMeanX << endl; *fLog << " - Meany [mm] = " << fMeanY << endl; *fLog << " - atg(y/x) [deg] = " << atg << endl; } // -------------------------------------------------------------------------- // // Print the hillas Parameters to *fLog, depending on the geometry in // units of deg // void MHillas::Print(const MGeomCam &geom) const { const Double_t atg = atan2(fMeanX, fMeanY)*kRad2Deg; *fLog << all; *fLog << GetDescriptor() << endl; *fLog << " - Length [deg] = " << fLength*geom.GetConvMm2Deg() << endl; *fLog << " - Width [deg] = " << fWidth *geom.GetConvMm2Deg() << endl; *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl; *fLog << " - Size [phe] = " << fSize << endl; *fLog << " - Meanx [deg] = " << fMeanX *geom.GetConvMm2Deg() << endl; *fLog << " - Meany [deg] = " << fMeanY *geom.GetConvMm2Deg() << endl; *fLog << " - atg(y/x) [deg] = " << atg << endl; } // -------------------------------------------------------------------------- // // Paint the ellipse corresponding to the parameters // void MHillas::Paint(Option_t *opt) { static const Float_t sOffsetW=40.0; static const Float_t sOffsetL=266.0; if (fLength<0 || fWidth<0) return; const Float_t l = fLength+sOffsetL; const Float_t w = fWidth +sOffsetW; TLine line; line.SetLineWidth(1); line.SetLineColor(kRed); // Length line line.PaintLine(-l*fCosDelta+fMeanX, -l*fSinDelta+fMeanY, +l*fCosDelta+fMeanX, +l*fSinDelta+fMeanY); // Width line line.PaintLine(+w*fSinDelta+fMeanX, -w*fCosDelta+fMeanY, -w*fSinDelta+fMeanX, +w*fCosDelta+fMeanY); // Rough estimate for disp const Double_t p = 362*(1-fWidth/fLength); // Disp line line.PaintLine( p*fCosDelta+20*fSinDelta+fMeanX, p*fSinDelta-20*fCosDelta+fMeanY, p*fCosDelta-20*fSinDelta+fMeanX, p*fSinDelta+20*fCosDelta+fMeanY); line.PaintLine(-p*fCosDelta+20*fSinDelta+fMeanX, -p*fSinDelta-20*fCosDelta+fMeanY, -p*fCosDelta-20*fSinDelta+fMeanX, -p*fSinDelta+20*fCosDelta+fMeanY); TEllipse e(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta*kRad2Deg+180); e.SetLineWidth(2); e.SetLineColor(kGreen); e.Paint(); } // -------------------------------------------------------------------------- // // 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. // Returns: // 0 no error // 1 number of pixels in event == 0 (special MC events) // 2 size==0 // 3 number of used pixel < 3 // 4 CorrXY == 0 // Int_t MHillas::Calc(const MGeomCam &geom, const MSignalCam &evt, Int_t island) { const UInt_t numpix = evt.GetNumPixels(); // // sanity check 1 (special MC events) // if (numpix==0) return 1; // // 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; fMeanY = 0; fSize = 0; UInt_t numused = 0; for (UInt_t i=0; i=0 && pix.GetIdxIsland()!=island) continue; const Float_t nphot = pix.GetNumPhotons(); // [#phot] if (nphot<0) continue; const MGeom &gpix = geom[i]; fSize += nphot; // [counter] fMeanX += nphot * gpix.GetX(); // [mm] fMeanY += nphot * gpix.GetY(); // [mm] numused++; } // // sanity check 2 // if (fSize==0) return 2; fMeanX /= fSize; // [mm] fMeanY /= fSize; // [mm] // // sanity check 3 // if (numused<3) return 3; // // calculate 2nd moments // --------------------- // // To make sure that the more complicated sum is evaluated as // accurate as possible (because it is needed for more // complicated calculations (see below) we calculate it using // double prcision. // Double_t corrxx=0; // [m^2] Double_t corrxy=0; // [m^2] Double_t corryy=0; // [m^2] for (UInt_t i=0; i=0 && pix.GetIdxIsland()!=island) continue; const Float_t nphot = pix.GetNumPhotons(); // [#phot] if (nphot<0) continue; const MGeom &gpix = geom[i]; const Float_t dx = gpix.GetX() - fMeanX; // [mm] const Float_t dy = gpix.GetY() - fMeanY; // [mm] corrxx += nphot * dx*dx; // [mm^2] corrxy += nphot * dx*dy; // [mm^2] corryy += nphot * dy*dy; // [mm^2] } // // calculate the basic Hillas parameters: orientation and size of axes // ------------------------------------------------------------------- // // 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 // // Rem: I tested replacing sqrt() by hypot() but they exactly // consume the same amount of time // const Double_t d0 = corryy - corrxx; const Double_t d1 = corrxy*2; const Double_t d2 = TMath::Sqrt(d0*d0 + d1*d1) + d0; // [0 const Double_t tand = d2==0 ? 0 : d2 / d1; // Force 0 also if d1==0 const Double_t tand2 = tand*tand; const Double_t s2 = tand2+1; const Double_t s = TMath::Sqrt(s2); // Set default for the case in which the image is symmetric on the y-axis fDelta = TMath::Pi()/2; fCosDelta = 0; fSinDelta = 1; Double_t axis1 = corryy; Double_t axis2 = corrxx; // This are all cases in which the image is not symmetric on the y-axis if (d1!=0 || d2==0) { fDelta = TMath::ATan(tand); fCosDelta = 1.0 /s; // need these in derived classes fSinDelta = tand/s; // like MHillasExt axis1 = (tand2*corryy + d2 + corrxx)/s2; axis2 = (tand2*corrxx - d2 + corryy)/s2; } // // fLength^2 is the second moment along the major axis of the distribution // fWidth^2 is the second moment along the minor axis of the distribution // // From the algorithm we get: fWidth <= fLength is always true // // very small numbers can get negative by rounding // fLength = axis1<0 ? 0 : TMath::Sqrt(axis1/fSize); // [mm] fWidth = axis2<0 ? 0 : TMath::Sqrt(axis2/fSize); // [mm] SetReadyToSave(); return 0; } // -------------------------------------------------------------------------- // // This function is ment for special usage, please never try to set // values via this function // void MHillas::Set(const TArrayF &arr) { if (arr.GetSize() != 6) return; fLength = arr.At(0); // [mm] major axis of ellipse fWidth = arr.At(1); // [mm] minor axis of ellipse fDelta = arr.At(2); // [rad] angle of major axis with x-axis fSize = arr.At(3); // [#CerPhot] sum of content of all pixels (number of Cherenkov photons) fMeanX = arr.At(4); // [mm] x-coordinate of center of ellipse fMeanY = arr.At(5); // [mm] y-coordinate of center of ellipse }