Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 7176) +++ /trunk/MagicSoft/Mars/Changelog (revision 7177) @@ -43,4 +43,13 @@ * mmuon/MMuonCalibPar.cc, mmuon/MMuonSearchPar.cc: - fixed Print output + + * mimage/MHillas.cc: + - implemented plotting a estimate of disp + - fixed some units in output + - fixed the special case in which the shower is symmetric to + the axis + + * mimage/MHillasCalc.cc: + - skipped printing the removed case 4 from MHillas::Calc Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 7176) +++ /trunk/MagicSoft/Mars/NEWS (revision 7177) @@ -14,4 +14,6 @@ - general: MH::SetPalette offers a lot of new palettes + + - general: MHillas - the case of CorrXY==0 is now handled properly - showplot: got a new option to start a root interpreter, too Index: /trunk/MagicSoft/Mars/mimage/MHillas.cc =================================================================== --- /trunk/MagicSoft/Mars/mimage/MHillas.cc (revision 7176) +++ /trunk/MagicSoft/Mars/mimage/MHillas.cc (revision 7177) @@ -159,5 +159,5 @@ *fLog << " - Width [mm] = " << fWidth << endl; *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl; - *fLog << " - Size [1] = " << fSize << " #CherPhot" << endl; + *fLog << " - Size [phe] = " << fSize << endl; *fLog << " - Meanx [mm] = " << fMeanX << endl; *fLog << " - Meany [mm] = " << fMeanY << endl; @@ -182,5 +182,5 @@ *fLog << " - Width [deg] = " << fWidth *geom.GetConvMm2Deg() << endl; *fLog << " - Delta [deg] = " << fDelta*kRad2Deg << endl; - *fLog << " - Size [1] = " << fSize << " #CherPhot" << endl; + *fLog << " - Size [phe] = " << fSize << endl; *fLog << " - Meanx [deg] = " << fMeanX *geom.GetConvMm2Deg() << endl; *fLog << " - Meany [deg] = " << fMeanY *geom.GetConvMm2Deg() << endl; @@ -215,4 +215,13 @@ 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); @@ -259,9 +268,5 @@ UInt_t numused = 0; -/* - MSignalPix *pix = 0; - - TIter Next(evt); - while ((pix=(MSignalPix*)Next()))*/ + for (UInt_t i=0; iGetPixId()*/]; + const MGeomPix &gpix = geom[i]; const Float_t nphot = pix.GetNumPhotons(); @@ -312,7 +317,4 @@ Double_t corryy=0; // [m^2] - //Next.Reset(); - //while ((pix=(MSignalPix*)Next())) - //{ for (UInt_t i=0; iGetPixId()*/]; + const MGeomPix &gpix = geom[i]; const Float_t dx = gpix.GetX() - fMeanX; // [mm] @@ -337,15 +339,4 @@ // - // If corrxy=0 (which should happen not very often, because fSize>0) - // we cannot calculate Length and Width. - // In reallity it is almost impossible to have a distribution of - // cerenkov photons in the used pixels which is exactly symmetric - // along one of the axis. - // It seems to be less than 0.1% of all events. - // - if (corrxy==0) - return 4; - - // // calculate the basic Hillas parameters: orientation and size of axes // ------------------------------------------------------------------- @@ -360,22 +351,35 @@ const Double_t d0 = corryy - corrxx; const Double_t d1 = corrxy*2; - const Double_t d2 = d0 + TMath::Sqrt(d0*d0 + d1*d1); - const Double_t tand = d2 / d1; + 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; - fDelta = TMath::ATan(tand); - - const Double_t s2 = tand2+1; - const Double_t s = TMath::Sqrt(s2); - - fCosDelta = 1.0/s; // need these in derived classes - fSinDelta = tand/s; // like MHillasExt - - const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/fSize; - const Double_t axis2 = (tand2*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 + 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 @@ -383,6 +387,6 @@ // very small numbers can get negative by rounding // - fLength = axis1<0 ? 0 : TMath::Sqrt(axis1); // [mm] - fWidth = axis2<0 ? 0 : TMath::Sqrt(axis2); // [mm] + fLength = axis1<0 ? 0 : TMath::Sqrt(axis1/fSize); // [mm] + fWidth = axis2<0 ? 0 : TMath::Sqrt(axis2/fSize); // [mm] SetReadyToSave(); Index: /trunk/MagicSoft/Mars/mimage/MHillasCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mimage/MHillasCalc.cc (revision 7176) +++ /trunk/MagicSoft/Mars/mimage/MHillasCalc.cc (revision 7177) @@ -376,5 +376,5 @@ PrintSkipped(fErrors[2], "Calculated Size == 0 (after cleaning)"); PrintSkipped(fErrors[3], "Number of used pixels < 3"); - PrintSkipped(fErrors[4], "CorrXY==0"); + //PrintSkipped(fErrors[4], "CorrXY==0"); } if (TestFlag(kCalcHillasSrc))