Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 7190) +++ trunk/MagicSoft/Mars/NEWS (revision 7191) @@ -99,4 +99,7 @@ - ganymed: optimized palettes for MHDisp + + - ganymed: MPedestalY3 now uses the same extraction range as the data + as calculated by MPedestalY2 - ganymed: Fixed an important bug when using ThetaSq as fit function: Index: trunk/MagicSoft/Mars/mcalib/MCalibrationRelTimeCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationRelTimeCalc.cc (revision 7190) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationRelTimeCalc.cc (revision 7191) @@ -472,81 +472,56 @@ void MCalibrationRelTimeCalc::FinalizeUnsuitablePixels() { - - *fLog << inf << endl; - *fLog << GetDescriptor() << ": Rel. Times Calibration status:" << endl; - *fLog << dec << setfill(' '); - - MCalibrationRelTimeCam *relcam = fIntensCam - ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam; - MBadPixelsCam *badcam = fIntensBad - ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels; - - const Int_t nareas = fGeom->GetNumAreas(); - - TArrayI counts(nareas); - - for (Int_t i=0; iGetSize(); i++) - { - MBadPixelsPix &bad = (*badcam)[i]; - if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun)) + *fLog << inf << endl; + *fLog << GetDescriptor() << ": Rel. Times Calibration status:" << endl; + *fLog << dec; + + MCalibrationRelTimeCam *relcam = fIntensCam ? (MCalibrationRelTimeCam*)fIntensCam->GetCam() : fCam; + const MBadPixelsCam *badcam = fIntensBad ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels; + + const Int_t nareas = fGeom->GetNumAreas(); + + TArrayI unsuit(nareas); + TArrayI unrel(nareas); + + for (int aidx=0; aidxGetNumUnsuitable(MBadPixelsPix::kUnsuitableRun, fGeom, aidx); + unrel[aidx] += badcam->GetNumUnsuitable(MBadPixelsPix::kUnreliableRun, fGeom, aidx); + relcam->SetNumUnsuitable(unsuit[aidx], aidx); + relcam->SetNumUnreliable(unrel[aidx], aidx); + } + + if (fGeom->InheritsFrom("MGeomCamMagic")) + { + *fLog << " Uncalibrated Pixels: Inner: " + << Form("%3i",unsuit[0]) << " Outer: " << Form("%3i",unsuit[1]) << endl; + *fLog << " Unreliable Pixels: Inner: " + << Form("%3i",unrel[0]) << " Outer: " << Form("%3i",unrel[1]) << endl; + } +} + +// ----------------------------------------------------------------------------------------------- +// +// Print out statistics about BadPixels of type UncalibratedType_t +// +void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const +{ + const MBadPixelsCam *badcam = fIntensBad ? (MBadPixelsCam*)fIntensBad->GetCam() : fBadPixels; + + UInt_t countinner = 0; + UInt_t countouter = 0; + for (Int_t i=0; iGetSize(); i++) + { + if ((*badcam)[i].IsUncalibrated(typ)) { - const Int_t aidx = (*fGeom)[i].GetAidx(); - counts[aidx]++; + if (fGeom->GetPixRatio(i) == 1.) + countinner++; + else + countouter++; } } - for (Int_t aidx=0; aidxSetNumUnsuitable(counts[aidx], aidx); - - if (fGeom->InheritsFrom("MGeomCamMagic")) - *fLog << " " << setw(7) << "Uncalibrated Pixels: Inner: " - << Form("%3i",counts[0]) << " Outer: " << Form("%3i",counts[1]) << endl; - - counts.Reset(); - - for (Int_t i=0; iGetSize(); i++) - { - MBadPixelsPix &bad = (*badcam)[i]; - if (bad.IsUnsuitable(MBadPixelsPix::kUnreliableRun)) - { - const Int_t aidx = (*fGeom)[i].GetAidx(); - counts[aidx]++; - } - } - - for (Int_t aidx=0; aidxSetNumUnreliable(counts[aidx], aidx); - - *fLog << " " << setw(7) << "Unreliable Pixels: Inner: " - << Form("%3i",counts[0]) << " Outer: " << Form("%3i",counts[1]) << endl; - -} - -// ----------------------------------------------------------------------------------------------- -// -// Print out statistics about BadPixels of type UncalibratedType_t -// -void MCalibrationRelTimeCalc::PrintUncalibrated(MBadPixelsPix::UncalibratedType_t typ, const char *text) const -{ - - MBadPixelsCam *badcam = fIntensBad - ? (MBadPixelsCam*) fIntensBad->GetCam() : fBadPixels; - - UInt_t countinner = 0; - UInt_t countouter = 0; - for (Int_t i=0; iGetSize(); i++) - { - MBadPixelsPix &bad = (*badcam)[i]; - if (bad.IsUncalibrated(typ)) - { - if (fGeom->GetPixRatio(i) == 1.) - countinner++; - else - countouter++; - } - } - - *fLog << " " << setw(7) << text - << Form("%s%3i%s%3i","Inner: ",countinner," Outer: ",countouter) << endl; + *fLog << " " << text << "Inner: " << Form("%3i",countinner); + *fLog << " Outer: " << Form("%3i", countouter) << endl; } Index: trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 7190) +++ trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 7191) @@ -27,4 +27,59 @@ // MFMagicCuts // +// Predefinitions: +// --------------- +// +// m3long = MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*fMM2Deg +// antim3long = MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMM2Deg +// +// Coefficients/Cuts: +// ------------------ +// +// c[0], c[5], c[6], c[7]: +// 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]) +// thetasq = disp^2 + dist^2 - 2*disp*dist*alpha +// antithetasq = antidisp^2 + dist^2 - 2*antidisp*dist*alpha +// +// c[1]: +// thetasq < c[1]*c[1] +// +// c[2], c[3], c[4]: +// A = c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3])) +// area < A +// +// +// Options: +// -------- +// +// HadronnessCut: +// - none/nocut: No area cut +// - area: Area cut +// - all: same as area +// +// ThetaCut: +// - none/nocut: no theta cut +// - on: cut in theta only +// - off: anti-cut in anti-theta only +// - wobble: same as on|off (both) +// +// +// Input Container: +// ------ +// MGeomCam +// MHillas +// MHillasExt +// MNewImagePar +// MHillasSrc +// [MHillasSrcAnti [MHillasSrc]] +// +// Output: +// ------- +// ThetaSquared [MParameterD] // stores thetasq +// Disp [MParameterD] // stores -p*sign(MHillasSrc.fCosDeltaAlpha) +// ThetaSquaredCut [MParameterD] // stores c[1]*c[1] +// ///////////////////////////////////////////////////////////////////////////// #include "MFMagicCuts.h" @@ -36,4 +91,6 @@ #include "MLog.h" #include "MLogManip.h" + +#include "MMath.h" #include "MParList.h" @@ -42,6 +99,6 @@ #include "MHillasSrc.h" #include "MHillasExt.h" +#include "MNewImagePar.h" #include "MParameters.h" -#include "MPointingPos.h" ClassImp(MFMagicCuts); @@ -54,6 +111,6 @@ // MFMagicCuts::MFMagicCuts(const char *name, const char *title) - : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fThetaSq(0), - fDisp(0), fPointing(0), fMatrix(0), fVariables(30), fThetaCut(kNone), + : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fNewImgPar(0), + fThetaSq(0), fDisp(0), fMatrix(0), fVariables(30), fThetaCut(kNone), fHadronnessCut(kAll) { @@ -68,7 +125,7 @@ AddToBranchList("MHillasSrc.fAlpha"); AddToBranchList("MHillasSrcAnti.fAlpha"); - AddToBranchList("MPointingPos.fZd"); AddToBranchList("MHillasExt.fMaxDist"); AddToBranchList("MHillasExt.fM3Long"); + AddToBranchList("MNewImagePar.fLeakage1"); fVariables[0] = 1.42547; // Xi @@ -110,20 +167,4 @@ return kFALSE; thetacut->SetVal(GetThetaSqCut()); - //thetacut->SetReadyToSave(); - - MParameterD *m3lcut = (MParameterD*)pList->FindCreateObj("MParameterD", "M3lCut"); - if (!m3lcut) - return kFALSE; - m3lcut->SetVal(fVariables[5]); - - MParameterD *dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXi"); - if (!dispxi) - return kFALSE; - dispxi->SetVal(fVariables[0]); - - dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXiTheta"); - if (!dispxi) - return kFALSE; - dispxi->SetVal(fVariables[6]); Print(); @@ -148,8 +189,8 @@ } - fPointing = (MPointingPos*)pList->FindObject("MPointingPos"); - if (!fPointing) - { - *fLog << err << "MPointingPos not found... aborting." << endl; + fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar"); + if (!fNewImgPar) + { + *fLog << err << "MNewImagePar not found... aborting." << endl; return kFALSE; } @@ -200,23 +241,20 @@ fMatrix = mat; - fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)"); - fMap[kEArea] = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg"); - fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg"); - fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength"); - fMap[kEZd] = fMatrix->AddColumn("MPointingPos.GetZdRad"); - - fMap[kEAlpha] = fMatrix->AddColumn("MHillasSrc.fAlpha"); - fMap[kEDist] = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg"); - //fMap[kDCA] = fMatrix->AddColumn("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg"); + fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)"); + fMap[kEArea] = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg"); + fMap[kEM3Trans] = fMatrix->AddColumn("abs(MHillasExt.fM3Trans)*MGeomCam.fConvMm2Deg"); + fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg"); + fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)"); + fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength"); + + fMap[kEAlpha] = fMatrix->AddColumn("MHillasSrc.fAlpha"); + fMap[kEDist] = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg"); + fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)"); + if (fThetaCut&kOff) { fMap[kEAlphaAnti] = fMatrix->AddColumn("MHillasSrcAnti.fAlpha"); fMap[kEDistAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg"); - fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg"); - //fMap[kDCAAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDCA*MGeomCam.fConvMm2Deg"); - } - - //fMap[kLength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg"); - //fMap[kWidth] = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg"); + } } @@ -229,13 +267,10 @@ // If par!=NULL p is stored in this MParameterD // -Double_t MFMagicCuts::GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par) const +Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const { // wl = width/l [1] // d = dist [deg] // a = alpha [deg] - const Double_t p = c*(1-wl); - if (par) - par->SetVal(p); - return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad()); + return d*d + p*p - 2*d*p*TMath::Cos(a*TMath::DegToRad()); } @@ -256,14 +291,16 @@ { // Get some variables - const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength(); - const Double_t lgsize = fMatrix ? GetVal(kESize) : TMath::Log10(fHil->GetSize()); - const Double_t zd = fMatrix ? GetVal(kEZd) : fPointing->GetZdRad(); + const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength(); + const Double_t lgsize = fMatrix ? GetVal(kESize) : TMath::Log10(fHil->GetSize()); + //const Double_t zd = fMatrix ? GetVal(kEZd) : fPointing->GetZdRad(); + const Double_t leak = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1); // For simplicity const Double_t *c = fVariables.GetArray(); - // Value for Theta cut - const Double_t cut = GetThetaSqCut(); //c[1]*c[1]; - const Double_t xi = (c[0]+c[8]*(lgsize-c[7])*(lgsize-c[7])) - c[6]*zd*zd; + // Value for Theta cut (Disp parametrization) + const Double_t cut = GetThetaSqCut(); + const Double_t xi = c[0]+c[6]*pow(leak, c[7]); + const Double_t disp = xi*(1-wdivl); // Default if we return before the end @@ -272,11 +309,20 @@ // ------------------- Most efficient cut ----------------------- // ---------------------- Theta cut ON -------------------------- - const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg; - const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha(); - - const Double_t thetasq = GetThetaSq(xi, wdivl, dist, alpha, fDisp); - + const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg; + const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha(); + 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()); + + // Align disp along source direction depending on third moment + const Double_t p = TMath::Sign(disp, m3long-c[5]); + + // Align sign of disp along shower axis like m3long + fDisp->SetVal(-p*sign); + + // Calculate corresponding Theta^2 + const Double_t thetasq = GetThetaSq(p, wdivl, dist, alpha); fThetaSq->SetVal(thetasq); + // Perform theta cut if (fThetaCut&kOn && thetasq >= cut) return kTRUE; @@ -290,23 +336,29 @@ if (area>=A) return kTRUE; - } - if (fHadronnessCut&kM3Long) - { - const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha()); - if (m3long<=c[5]) - return kTRUE; + + //const Double_t m3t = fMatrix ? GetVal(kEM3Trans) : TMath::Abs(fHilExt->GetM3Trans())*fMm2Deg; + //if (m3t>c[8]) + // return kTRUE; } // ------------------- Least efficient cut ----------------------- // ---------------------- Theta cut OFF -------------------------- + if (fThetaCut&kOff) { const Double_t distanti = fMatrix ? GetVal(kEDistAnti) : fHilAnti->GetDist()*fMm2Deg; const Double_t alphaanti = fMatrix ? GetVal(kEAlphaAnti) : fHilAnti->GetAlpha(); - const Double_t thetasqanti = GetThetaSq(xi, wdivl, distanti, alphaanti); const Double_t m3lanti = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha()); - if (thetasqantic[5])) + const Double_t panti = TMath::Sign(disp, m3lanti-c[5]); + const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti); + + if (thetasqanti