/* ======================================================================== *\ ! ! * ! * 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): Rudolf Bock 10/2001 ! Author(s): Wolfgang Wittek 06/2002 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MHillasExt // // Storage Container for extended image parameters // // extended image parameters // // // Version 1: // ---------- // fConc ratio of sum of two highest pixels over fSize // fConc1 ratio of highest pixel over fSize // fAsym distance from highest pixel to center, projected onto major axis // fM3Long third moment along major axis // fM3Trans third moment along minor axis // // Version 2: // ---------- // fConc removed // fConc1 removed // // Version 3: // ---------- // fMaxDist added. Distance between center and most distant used pixel // // // WARNING: Before you can use fAsym, fM3Long and fM3Trans you must // multiply by the sign of MHillasSrc::fCosDeltaAlpha // //////////////////////////////////////////////////////////////////////////// /* // fAsymna d/(d na) of ( sum(x*q^na)/sum(q^na), sum(y*q^na)/sum(q^na) ) // projected onto the major axis // fAsym0 (F-B)/(F+B) along the major axis */ #include "MHillasExt.h" #include #include #include #include #include "MGeomPix.h" #include "MGeomCam.h" #include "MSignalPix.h" #include "MSignalCam.h" #include "MLog.h" #include "MLogManip.h" #include "MHillas.h" ClassImp(MHillasExt); using namespace std; // ------------------------------------------------------------------------- // // Default constructor. // MHillasExt::MHillasExt(const char *name, const char *title) { fName = name ? name : "MHillasExt"; fTitle = title ? title : "Storage container for extended parameter set of one event"; Reset(); } // ------------------------------------------------------------------------- // // Reset values to: // fAsym = 0; // fM3Long = 0; // fM3Trans = 0; // fMaxDist = 0; // void MHillasExt::Reset() { fAsym = 0; fM3Long = 0; fM3Trans = 0; fMaxDist = 0; } // ------------------------------------------------------------------------- // // calculation of additional parameters based on the camera geometry // and the cerenkov photon event // Int_t MHillasExt::Calc(const MGeomCam &geom, const MSignalCam &evt, const MHillas &hil, Int_t island) { // // calculate the additional image parameters // -------------------------------------------- // // loop to get third moments along ellipse axes and two max pixels // // For the moments double precision is used to make sure, that // the complex matrix multiplication and sum is evaluated correctly. // Double_t m3x = 0; Double_t m3y = 0; Int_t maxpixid = 0; Float_t maxpix = 0; Float_t maxdist = 0; // MSignalPix *pix = 0; // TIter Next(evt); // while ((pix=(MSignalPix*)Next())) const UInt_t npix = evt.GetNumPixels(); for (UInt_t i=0; i=0 && pix.GetIdxIsland()!=island) continue; //const Int_t pixid = pix->GetPixId(); const MGeomPix &gpix = geom[i/*pixid*/]; const Double_t dx = gpix.GetX() - hil.GetMeanX(); // [mm] const Double_t dy = gpix.GetY() - hil.GetMeanY(); // [mm] Double_t nphot = pix.GetNumPhotons(); // [1] const Double_t dzx = hil.GetCosDelta()*dx + hil.GetSinDelta()*dy; // [mm] const Double_t dzy = -hil.GetSinDelta()*dx + hil.GetCosDelta()*dy; // [mm] const Double_t dist = dx*dx+dy*dy; if (dist>TMath::Abs(maxdist)) maxdist = dzx<0 ? -dist : dist; // [mm^2] m3x += nphot * dzx*dzx*dzx; // [mm^3] m3y += nphot * dzy*dzy*dzy; // [mm^3] // // Now we are working on absolute values of nphot, which // must take pixel size into account // nphot *= geom.GetPixRatio(i/*pixid*/); if (nphot>maxpix) { maxpix = nphot; // [1] maxpixid = i;//pixid; } } const MGeomPix &maxp = geom[maxpixid]; // // Asymmetry // fAsym = (hil.GetMeanX()-maxp.GetX())*hil.GetCosDelta() + (hil.GetMeanY()-maxp.GetY())*hil.GetSinDelta(); // [mm] // // Third moments along axes get normalized // m3x /= hil.GetSize(); m3y /= hil.GetSize(); fM3Long = m3x<0 ? -pow(-m3x, 1./3) : pow(m3x, 1./3); // [mm] fM3Trans = m3y<0 ? -pow(-m3y, 1./3) : pow(m3y, 1./3); // [mm] // // Distance between max signal and COG // const Double_t md = TMath::Sqrt(TMath::Abs(maxdist)); fMaxDist = maxdist<0 ? -md : md; // [mm] SetReadyToSave(); return 0; } // -------------------------------------------------------------------------- // // This function is ment for special usage, please never try to set // values via this function // void MHillasExt::Set(const TArrayF &arr) { if (arr.GetSize() != 4) return; fAsym = arr.At(0); fM3Long = arr.At(1); fM3Trans = arr.At(2); fMaxDist = arr.At(3); } // ------------------------------------------------------------------------- // // Print contents of MHillasExt to *fLog. // void MHillasExt::Print(Option_t *) const { *fLog << all; *fLog << GetDescriptor() << endl; *fLog << " - Asymmetry [mm] = " << fAsym << endl; *fLog << " - 3.Moment Long [mm] = " << fM3Long << endl; *fLog << " - 3.Moment Trans [mm] = " << fM3Trans << endl; *fLog << " - Max.Dist [mm] = " << fMaxDist << endl; } // ------------------------------------------------------------------------- // // Print contents of MHillasExt to *fLog, depending on the geometry in // units of deg. // void MHillasExt::Print(const MGeomCam &geom) const { *fLog << all; *fLog << GetDescriptor() << endl; *fLog << " - Asymmetry [deg] = " << fAsym *geom.GetConvMm2Deg() << endl; *fLog << " - 3.Moment Long [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl; *fLog << " - 3.Moment Trans [deg] = " << fM3Trans*geom.GetConvMm2Deg() << endl; *fLog << " - Max.Dist [deg] = " << fMaxDist*geom.GetConvMm2Deg() << endl; } // ------------------------------------------------------------------------- // // Paint the 3rd moment on top of the shower. Therefor "MHillas" is needed. // it is searched via gPad->FindObject. If opt is not IsNull an object with // the name opt is searched. // void MHillasExt::Paint(Option_t *opt) { const TString name(opt); const MHillas *hil = dynamic_cast(gPad->FindObject(name.IsNull() ? "MHillas" : name.Data())); if (!hil) return; TVector2 v(fM3Long, fM3Trans); v = v.Rotate(hil->GetDelta()+TMath::Pi()); v += hil->GetMean(); TMarker m; m.SetMarkerColor(15); m.SetMarkerStyle(kFullDotLarge); m.PaintMarker(v.X(), v.Y()); }