Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 8618) +++ /trunk/MagicSoft/Mars/Changelog (revision 8619) @@ -20,4 +20,19 @@ + 2007/06/29 Thomas Bretz + + * mimage/MHillasExt.[h,cc]: + - added new timing parameters fSlopeTrans and fSlopeLong + - removed never used parameter fMaxDist + - increased accordingly the class version number by 1 + - replaced the old calculation of the third moments by a + numerically more accurate one, by avoiding to calculate + to many differences too often. + + * mfilter/MFMagicCuts.[h,cc]: + - added the usage of the new SlopeLong parameter for ghostbusting + + + 2007/06/28 Thomas Bretz @@ -48,4 +63,24 @@ * mjobs/MJCalibrateSignal.cc: - changed axis label of PulsePos plot (now in nanosec) + + * mpointing/MSrcPosCalc.[h,cc]: + - allow to set a tasklist as callback to now which n-th + pass of the same task list it is + + * mbase/MTaskList.[h,cc]: + - added some code to allow the execution of one task list more + than once. This is for example necessary to process three + different off-source regions. + + * mjobs/MJCut.[h,cc]: + - use the new feature in MTaskList to setup a tasklist + processing the off-source calculation tasklist more than once + - added a new data meber fNumOffSourcePos + - added a new resource option NumOffSourcePositions + - added a new CutQ before Cut0 which takes place before all + source posisiton dependant stuff + + * ganymed_onoff.rc: + - renamed Cut0 to CutQ Index: /trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 8618) +++ /trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 8619) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 03/2005 ! -! Copyright: MAGIC Software Development, 2000-2005 +! Copyright: MAGIC Software Development, 2000-2007 ! ! @@ -38,8 +38,10 @@ // dist = MHillasSrc.fDist*fMm2Deg // m3long = MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*fMm2Deg +// slope = MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)*fMm2Deg // // antialpha = MHillasSrcAnti.fAlpha // antidist = MHillasSrcAnti.fDist*fMm2Deg // antim3long = MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMm2Deg +// antislope = MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMm2Deg // // had = Hadronness.fVal @@ -48,9 +50,13 @@ // ------------------ // -// c[0], c[5], c[6], c[7]: +// c[0], c[5], c[6], c[7], c[8], c[9]: // xi = c[0]+c[6]*pow(leakage1, c[7]); // p = xi*(1-width/length); -// disp = TMath::Sign(disp, m3long-c[5]) -// antidisp = TMath::Sign(disp, antim3long-c[5]) +// sign1 = m3long-c[5] +// sign2 = (dist-c[9])*c[8]-slope +// disp = sign1<0 ||sign2<0 ? -disp : disp +// antisign1 = antim3long-c[5] +// antisign2 = (antidist-c[9])*c[8]-antislope +// antidisp = antisign1<0 || antisign2<0 ? -antidisp : antidisp // thetasq = disp^2 + dist^2 - 2*disp*dist*alpha // antithetasq = antidisp^2 + antidist^2 - 2*antidisp*antidist*antialpha @@ -65,7 +71,7 @@ // // HadronnessCut: -// c[8], c[9]: -// had <= c[8] -// 10^lgsize >= c[9] +// c[10], c[11]: +// had <= c[10] +// 10^lgsize >= c[11] // // @@ -150,4 +156,5 @@ AddToBranchList("MHillasSrcAnti.fCosDeltaAlpha"); AddToBranchList("MHillasExt.fM3Long"); + AddToBranchList("MHillasExt.fSlopeLong"); AddToBranchList("MNewImagePar.fLeakage1"); AddToBranchList("Hadronness.fVal"); @@ -277,4 +284,7 @@ fMatrix = mat; +// fMap[kELength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg"); +// fMap[kEWidth] = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg"); + fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength"); fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)"); @@ -289,4 +299,6 @@ fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)"); + fMap[kESlope] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg"); + if (fThetaCut&kOff) { @@ -294,4 +306,5 @@ fMap[kEDistAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg"); fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg"); + fMap[kESlopeAnti] = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg"); } @@ -351,9 +364,15 @@ const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg; const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha(); + const Double_t sign = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha()); const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha()); - const Double_t sign = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha()); + const Double_t slope = fMatrix ? GetVal(kESlope) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha()); + + // Do Ghostbusting with slope and m3l + const Double_t sign1 = (dist-c[9])*c[8]-slope; + const Double_t sign2 = m3long-c[5]; + const Double_t p = sign1<0 || sign2<0 ? -disp : disp; // Align disp along source direction depending on third moment - const Double_t p = TMath::Sign(disp, m3long-c[5]); + //const Double_t p = TMath::Sign(disp, m3long-c[5]); // Align sign of disp along shower axis like m3long @@ -384,8 +403,8 @@ { const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal(); - if (had>c[8]) + if (had>c[10]) return kTRUE; - if (TMath::Power(10, lgsize)GetAlpha(); const Double_t m3lanti = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha()); - - const Double_t panti = TMath::Sign(disp, m3lanti-c[5]); + const Double_t slopeanti = fMatrix ? GetVal(kESlopeAnti) : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha()); + + // Do Ghostbusting with slope and m3l + const Double_t sign3 = (distanti-c[9])*c[8]-slopeanti; + const Double_t sign4 = m3lanti-c[5]; + const Double_t panti = sign3<0 || sign4<0 ? -disp : disp; + + // Align disp along source direction depending on third moment + //const Double_t panti = TMath::Sign(disp, m3lanti-c[5]); const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti); Index: /trunk/MagicSoft/Mars/mfilter/MFMagicCuts.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFMagicCuts.h (revision 8618) +++ /trunk/MagicSoft/Mars/mfilter/MFMagicCuts.h (revision 8619) @@ -43,5 +43,6 @@ enum { kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist, kEDistAnti, - kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESrcSign, kEHadronness, + kEM3Long, kEM3LongAnti, kEWdivL, kELeakage, kESrcSign, + kESlope, kESlopeAnti, kEHadronness, kLastElement }; Index: /trunk/MagicSoft/Mars/mimage/MHillasExt.cc =================================================================== --- /trunk/MagicSoft/Mars/mimage/MHillasExt.cc (revision 8618) +++ /trunk/MagicSoft/Mars/mimage/MHillasExt.cc (revision 8619) @@ -20,5 +20,5 @@ ! Author(s): Wolfgang Wittek 06/2002 ! -! Copyright: MAGIC Software Development, 2000-2004 +! Copyright: MAGIC Software Development, 2000-2007 ! ! @@ -51,6 +51,12 @@ // fMaxDist added. Distance between center and most distant used pixel // -// -// WARNING: Before you can use fAsym, fM3Long and fM3Trans you must +// Version 4: +// ---------- +// fMaxDist removed. +// fSlopeLong added +// fSlopeTrans added +// +// +// WARNING: Before you can use fAsym, fSlope* and fM3* you must // multiply by the sign of MHillasSrc::fCosDeltaAlpha // @@ -66,4 +72,5 @@ #include #include +#include #include @@ -97,17 +104,13 @@ // ------------------------------------------------------------------------- // -// Reset values to: -// fAsym = 0; -// fM3Long = 0; -// fM3Trans = 0; -// fMaxDist = 0; +// Reset all values to 0 // void MHillasExt::Reset() { - fAsym = 0; - fM3Long = 0; - fM3Trans = 0; - - fMaxDist = 0; + fAsym = 0; + fM3Long = 0; + fM3Trans = 0; + fSlopeLong = 0; + fSlopeTrans = 0; } @@ -128,14 +131,23 @@ // 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())) + + // Variables to caluclate time slope + // Formula: http://mo.mathematik.uni-stuttgart.de/inhalt/erlaeuterung/erlaeuterung300/ + UInt_t cnt = 0; + + Double_t sumx = 0; + Double_t sumy = 0; + Double_t sumt = 0; + Double_t sumxy = 0; + Double_t sumxt = 0; + Double_t sumyt = 0; + Double_t sumx2 = 0; + Double_t sumy2 = 0; + Double_t sumx3 = 0; + Double_t sumy3 = 0; + Double_t sumx2y = 0; + Double_t sumxy2 = 0; const UInt_t npix = evt.GetNumPixels(); @@ -149,21 +161,33 @@ 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] + const MGeomPix &gpix = geom[i]; + + const Double_t x = gpix.GetX(); + const Double_t y = gpix.GetY(); + const Double_t t = pix.GetArrivalTime(); + + // --- time slope ---- + sumx += x; + sumy += y; + sumt += t; + + sumx2 += x*x; + sumy2 += y*y; + sumxy += x*y; + sumxt += x*t; + sumyt += y*t; + + // --- 3rd moment --- + const Double_t dx = x - hil.GetMeanX(); // [mm] + const Double_t dy = y - hil.GetMeanY(); // [mm] + + Double_t nphot = pix.GetNumPhotons(); // [1] + + sumx3 += nphot * dx*dx*dx; + sumy3 += nphot * dy*dy*dy; + sumx2y += nphot * dx*dx*dy; + sumx2y += nphot * dx*dy*dy; + + cnt++; // @@ -171,35 +195,49 @@ // must take pixel size into account // - nphot *= geom.GetPixRatio(i/*pixid*/); - + nphot *= geom.GetPixRatio(i); + + // --- max pixel --- if (nphot>maxpix) { - maxpix = nphot; // [1] - maxpixid = i;//pixid; + maxpix = nphot; // [1] + maxpixid = i; } } + const Double_t c = hil.GetCosDelta(); + const Double_t s = hil.GetSinDelta(); + + // + // Time slope + // + const Double_t dxt = c*sumxt + s*sumyt; + const Double_t dyt = -s*sumxt + c*sumyt; + + const Double_t dx = c*sumx + s*sumy; + const Double_t dy = -s*sumx + c*sumy; + + const Double_t dx2 = c*c*sumx2 + 2*c*s*sumxy + s*s*sumy2; + const Double_t dy2 = s*s*sumx2 - 2*c*s*sumxy + c*c*sumy2; + + const Double_t detx = cnt*dx2 - dx*dx; + const Double_t dety = cnt*dy2 - dy*dy; + + fSlopeLong = detx==0 ? 0 : (cnt*dxt - sumt*dx)/detx; + fSlopeTrans = dety==0 ? 0 : (cnt*dyt - sumt*dy)/dety; + + // + // Third moments along axes get normalized + // + const Double_t m3l = c*c*c*sumx3 + 3*(s*c*c*sumx2y + c*s*s*sumxy2) + s*s*s*sumy3; + const Double_t m3t = c*c*c*sumy3 + 3*(s*s*c*sumx2y - s*c*c*sumxy2) - s*s*s*sumx3; + + fM3Long = MMath::Sqrt3(m3l/hil.GetSize()); // [mm] + fM3Trans = MMath::Sqrt3(m3t/hil.GetSize()); // [mm] + + // + // Asymmetry + // 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] + fAsym = (hil.GetMeanX()-maxp.GetX())*c + (hil.GetMeanY()-maxp.GetY())*s; // [mm] SetReadyToSave(); @@ -215,11 +253,12 @@ void MHillasExt::Set(const TArrayF &arr) { - if (arr.GetSize() != 4) + if (arr.GetSize() != 5) return; - fAsym = arr.At(0); - fM3Long = arr.At(1); - fM3Trans = arr.At(2); - fMaxDist = arr.At(3); + fAsym = arr.At(0); + fM3Long = arr.At(1); + fM3Trans = arr.At(2); + fSlopeLong = arr.At(3); + fSlopeTrans = arr.At(4); } @@ -228,12 +267,17 @@ // Print contents of MHillasExt to *fLog. // -void MHillasExt::Print(Option_t *) const -{ +void MHillasExt::Print(Option_t *o) const +{ + const Bool_t showtrans = TString(o).Contains("trans"); + *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; + if (showtrans) + *fLog << " - 3.Moment Trans [mm] = " << fM3Trans << endl; + *fLog << " - Slope Long [mm] = " << fSlopeLong << endl; + if (showtrans) + *fLog << " - Slope Trans [mm] = " << fSlopeTrans << endl; } @@ -249,6 +293,7 @@ *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; + *fLog << " - 3.Moment Trans [mm] = " << fM3Trans*geom.GetConvMm2Deg() << endl; + *fLog << " - Slope Long [mm] = " << fSlopeLong*geom.GetConvMm2Deg() << endl; + *fLog << " - Slope Trans [mm] = " << fSlopeTrans*geom.GetConvMm2Deg() << endl; } @@ -272,5 +317,5 @@ TMarker m; - m.SetMarkerColor(15); + m.SetMarkerColor(kRed); m.SetMarkerStyle(kFullDotLarge); m.PaintMarker(v.X(), v.Y()); Index: /trunk/MagicSoft/Mars/mimage/MHillasExt.h =================================================================== --- /trunk/MagicSoft/Mars/mimage/MHillasExt.h (revision 8618) +++ /trunk/MagicSoft/Mars/mimage/MHillasExt.h (revision 8619) @@ -20,5 +20,6 @@ Float_t fM3Trans; // [mm] 3rd moment (e-weighted) along minor axis - Float_t fMaxDist; // Distance between center and most distant used pixel + Float_t fSlopeLong; + Float_t fSlopeTrans; public: @@ -27,9 +28,9 @@ void Reset(); - Float_t GetAsym() const { return fAsym; } - Float_t GetM3Long() const { return fM3Long; } - Float_t GetM3Trans() const { return fM3Trans; } - - Float_t GetMaxDist() const { return fMaxDist; } + Float_t GetAsym() const { return fAsym; } + Float_t GetM3Long() const { return fM3Long; } + Float_t GetM3Trans() const { return fM3Trans; } + Float_t GetSlopeLong() const { return fSlopeLong; } + Float_t GetSlopeTrans() const { return fSlopeTrans; } Int_t Calc(const MGeomCam &geom, const MSignalCam &pix, @@ -43,5 +44,5 @@ void Set(const TArrayF &arr); - ClassDef(MHillasExt, 3) // Storage Container for extended Hillas Parameter + ClassDef(MHillasExt, 4) // Storage Container for extended Hillas Parameter }; #endif