Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7172) +++ trunk/MagicSoft/Mars/Changelog (revision 7173) @@ -21,4 +21,20 @@ -*-*- END OF LINE -*-*- + + 2005/06/28 Thomas Bretz + + * mfilter/MFMagicCuts.[h,cc]: + - first full implemtation... updates to come. + + * mhbase/MH.[h,cc]: + - added new member function to set several palettes + + * mhflux/MHDisp.[h,cc]: + - fixed z-axis dscription + - rotate filling off data by 180 deg. + - implemented subtracting off-data from on-data + - set different palettes + + 2005/06/27 Thomas Bretz Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 7172) +++ trunk/MagicSoft/Mars/NEWS (revision 7173) @@ -12,4 +12,6 @@ - general: MHCamera now displays the profiles in deg instead of mm + + - general: MH::SetPalette offers a lot of new palettes - callisto: MCalibrationHiLoCam can now be printed from its context @@ -27,4 +29,6 @@ argument + - ganymed: The first version of MFMagicCuts have been released + - ganymed: the Conc1 plot was incorrectly scaled in MHVsSize @@ -32,4 +36,6 @@ of bins for the signal region (NumBinsSignal) and the number of total bins (NumBInsTital) in the MHThetaSq histogram + + - ganymed: optimized palettes for MHDisp - sponde: the zenith angle distribution is now weighted instead of Index: trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 7172) +++ trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc (revision 7173) @@ -30,9 +30,23 @@ #include "MFMagicCuts.h" +#include + +#include "MHMatrix.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MHillasSrc.h" +#include "MHillasExt.h" +#include "MParameters.h" +#include "MPointingPos.h" + ClassImp(MFMagicCuts); using namespace std; - // -------------------------------------------------------------------------- // @@ -40,6 +54,427 @@ // 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), + fHadronnessCut(kAll) { fName = name ? name : "MFMagicCuts"; fTitle = title ? title : "Class to evaluate the MagicCuts"; -} + + fMatrix = NULL; + + AddToBranchList("MHillas.fWidth"); + AddToBranchList("MHillas.fLength"); + AddToBranchList("MHillas.fSize"); + AddToBranchList("MHillasSrc.fAlpha"); + AddToBranchList("MHillasSrcAnti.fAlpha"); + AddToBranchList("MPointingPos.fZd"); + AddToBranchList("MHillasExt.fMaxDist"); + AddToBranchList("MHillasExt.fM3Long"); + + fVariables[0] = 1.42547; // Xi + fVariables[1] = 0.233773; // Theta^2 + fVariables[2] = 0.2539460; // Area[0] + fVariables[3] = 5.2149800; // Area[1] + fVariables[4] = 0.1139130; // Area[2] + fVariables[5] = -0.0889; // M3long + fVariables[6] = 0.18; // Xi(theta) +} + +// -------------------------------------------------------------------------- +// +// The theta cut value GetThetaSqCut() is propageted to the parameter list +// as 'ThetaSquaredCut' as MParameterD so that it can be used somewhere else. +// +Int_t MFMagicCuts::PreProcess(MParList *pList) +{ + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared"); + if (!fThetaSq) + return kFALSE; + fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp"); + if (!fDisp) + return kFALSE; + + // propagate Theta cut to the parameter list + // so that it can be used somewhere else + MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut"); + if (!thetacut) + 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(); + + if (fMatrix) + return kTRUE; + + //----------------------------------------------------------- + + fHil = (MHillas*)pList->FindObject("MHillas"); + if (!fHil) + { + *fLog << err << "MHillas not found... aborting." << endl; + return kFALSE; + } + + fHilExt = (MHillasExt*)pList->FindObject("MHillasExt"); + if (!fHilExt) + { + *fLog << err << "MHillasExt not found... aborting." << endl; + return kFALSE; + } + + fPointing = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointing) + { + *fLog << err << "MPointingPos not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc"); + if (!fHilSrc) + { + *fLog << err << "MHillasSrc not found... aborting." << endl; + return kFALSE; + } + + if (fThetaCut&kOff) + { + fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc"); + if (!fHilAnti) + { + *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Returns the mapped value from the Matrix +// +Double_t MFMagicCuts::GetVal(Int_t i) const +{ + return (*fMatrix)[fMap[i]]; +} + +// -------------------------------------------------------------------------- +// +// You can use this function if you want to use a MHMatrix instead of the +// given containers. This function adds all necessary columns to the +// given matrix. Afterward you should fill the matrix with the corresponding +// data (eg from a file by using MHMatrix::Fill). If you now loop +// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process +// will take the values from the matrix instead of the containers. +// +void MFMagicCuts::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + 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"); + 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"); +} + +// -------------------------------------------------------------------------- +// +// Returns theta squared based on the formula: +// p := c*(1-width/length) +// return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad()); +// +// 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 +{ + // 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()); +} + +// -------------------------------------------------------------------------- +// +// Returns squared cut value in theta: fVariables[1]^2 +// +Double_t MFMagicCuts::GetThetaSqCut() const +{ + return fVariables[1]*fVariables[1]; +} + +// --------------------------------------------------------------------------- +// +// Evaluate dynamical magic-cuts +// +Int_t MFMagicCuts::Process() +{ + // 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(); + + // 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; + + // Default if we return before the end + fResult = kFALSE; + + // ------------------- 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); + + fThetaSq->SetVal(thetasq); + + if (fThetaCut&kOn && thetasq >= cut) + return kTRUE; + + // --------------------- Efficient cut ------------------------- + // ---------------------- Area/M3l cut -------------------------- + if (fHadronnessCut&kArea) + { + const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg; + const Double_t A = lgsize>c[3] ? c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3])); + 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; + } + + // ------------------- 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])) + return kTRUE; + } + + fResult = kTRUE; + + return kTRUE; +} + +void MFMagicCuts::SetVariables(const TArrayD &arr) +{ + fVariables = arr; +} + +TString MFMagicCuts::GetParam(Int_t i) const +{ + if (i>=fVariables.GetSize()) + return ""; + + return Form("Param[%2d]=%+f", i, fVariables[i]); +} + +void MFMagicCuts::Print(Option_t *o) const +{ + *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl; + + *fLog << "Using Theta cut: "; + switch (fThetaCut) + { + case kNone: + *fLog << "none" << endl; + break; + case kOn: + *fLog << "on" << endl; + break; + case kOff: + *fLog << "off" << endl; + break; + case kWobble: + *fLog << "on and off (wobble)" << endl; + break; + } + *fLog << "Using hadronness cut: "; + switch (fHadronnessCut) + { + case kNoCut: + *fLog << "none" << endl; + break; + case kArea: + *fLog << "area" << endl; + break; + case kM3Long: + *fLog << "m3long" << endl; + break; + case kAll: + *fLog << "all" << endl; + break; + } + *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl; + + const Int_t n = fVariables.GetSize(); + for (int i=0; i> fVariables[i]; + if (!fin) + { + *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl; + return kERROR; + } + } + return kTRUE; +} + +Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const +{ + ofstream fout(fname); + if (!fout) + { + *fLog << err << "Cannot open file " << fname << ": "; + *fLog << strerror(errno) << endl; + return kFALSE; + } + + for (int i=0; i +#endif + +class MParList; + +class MHillas; +class MHillasSrc; +class MHillasExt; +class MParameterD; +class MPointingPos; class MHMatrix; class MFMagicCuts : public MFilter { +public: + // Possible kind of theta cuts + enum ThetaCut_t { + kNone =BIT(0), + kOn =BIT(1), + kOff =BIT(2), + kWobble=kOn|kOff + }; + // Possible kind of hadronness cuts + enum HadronnessCut_t { + kNoCut =BIT(0), + kArea =BIT(1), + kM3Long=BIT(2), + kAll =kArea|kM3Long + }; + private: + // Elements for mapping. kLastElement must not be used and must be + // the last on in the list. It is used as counter for fMap. + enum { + kESize, kEAlpha, kEAlphaAnti, kEArea, kEDist, + kEM3Long, kEM3LongAnti, kEDistAnti, kEWdivL, kEZd, + kLastElement + }; - virtual Bool_t IsExpressionTrue() const { return kTRUE; } + MHillas *fHil; //! Pointer to MHillas container + MHillasSrc *fHilSrc; //! Pointer to MHillasSrc container + MHillasSrc *fHilAnti; //! Pointer to MHillasSrc container called MHillasSrcAnti + MHillasExt *fHilExt; //! Pointer to MHillasExt container + MParameterD *fThetaSq; //! Pointer to MParameterD container called ThetaSq + MParameterD *fDisp; //! Pointer to MParameterD container called Disp + MPointingPos *fPointing; //! Pointer to MPointingPos container + + Float_t fMm2Deg; //! Conversion factor from mm to deg, from MGeomCam + Bool_t fResult; //! Result of the filter evaluation + + Int_t fMap[kLastElement]; //! Mapping table for fast optimization + MHMatrix *fMatrix; //! Matrix thorugh which the mapped elements are accessed + + TArrayD fVariables; // Coefficients of cuts + + ThetaCut_t fThetaCut; // Which kind of theta cut should be evaluated + HadronnessCut_t fHadronnessCut; // Which kind of hadronness cut should be evaluated + + // MTask + Int_t PreProcess(MParList *pList); + Int_t Process(); + + // MFilter + Bool_t IsExpressionTrue() const { return fResult; } + + // MFMagicCuts + Double_t GetVal(Int_t i) const; + TString GetParam(Int_t i) const; + Double_t GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par=0) const; public: MFMagicCuts(const char *name=NULL, const char *title=NULL); - void InitMapping(MHMatrix *mat) { } - void StopMapping() { InitMapping(NULL); } + // Getter + Double_t GetThetaSqCut() const; - ClassDef(MFMagicCuts, 0) // A filter to evaluate the MagicCuts + // Setter + void SetThetaCut(ThetaCut_t c) { fThetaCut=c; } + void SetHadronnessCut(HadronnessCut_t c) { fHadronnessCut=c; } + + // MFMagicCuts + void InitMapping(MHMatrix *mat); + void StopMapping() { InitMapping(NULL); } + + Bool_t CoefficentsRead(const char *fname); + Bool_t CoefficentsWrite(const char *fname) const; + + // MParContainer + void SetVariables(const TArrayD &arr); + + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + + // TObject + void Print(Option_t *o="") const; + + ClassDef(MFMagicCuts, 1) // A filter to evaluate the MagicCuts }; Index: trunk/MagicSoft/Mars/mhbase/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.cc (revision 7172) +++ trunk/MagicSoft/Mars/mhbase/MH.cc (revision 7173) @@ -1204,5 +1204,5 @@ // -------------------------------------------------------------------------- // -// M.Gaug added this withouz Documentation +// // TH1I* MH::ProjectArray(const TArrayF &array, Int_t nbins, @@ -1251,5 +1251,5 @@ // -------------------------------------------------------------------------- // -// M.Gaug added this withouz Documentation +// // TH1I* MH::ProjectArray(const TArrayD &array, Int_t nbins, const char* name, const char* title) @@ -1297,5 +1297,5 @@ // -------------------------------------------------------------------------- // -// M.Gaug added this withouz Documentation +// // TH1I* MH::ProjectArray(const MArrayF &array, Int_t nbins, @@ -1307,5 +1307,5 @@ // -------------------------------------------------------------------------- // -// M.Gaug added this withouz Documentation +// // TH1I* MH::ProjectArray(const MArrayD &array, Int_t nbins, const char* name, const char* title) @@ -1324,2 +1324,141 @@ *fLog << "%) Evts skipped: " << str << endl; } + +// -------------------------------------------------------------------------- +// +// Calls gStyle->SetPalette. Allowed palettes are: +// pretty +// deepblue: darkblue -> lightblue +// lightblue: black -> blue -> white +// greyscale: black -> white +// glow1: black -> darkred -> orange -> yellow -> white +// glow2: +// glowsym: lightblue -> blue -> black -> darkred -> orange -> yellow -> white +// redish: darkred -> lightred +// bluish: darkblue -> lightblue +// small1: +// +// If the palette name contains 'inv' the order of the colors is inverted. +// +// The second argument determines the number of colors for the palette. +// The default is 50. 'pretty' always has 50 colors. +// +// (Remark: Drawing 3D object like TH2D with surf3 allows a maximum +// of 99 colors) +// +void MH::SetPalette(TString paletteName, Int_t ncol) +{ + Bool_t found=kFALSE; + + paletteName.ToLower(); + + const Bool_t inverse = paletteName.Contains("inv"); + + if (paletteName.Contains("pretty")) + { + gStyle->SetPalette(1, 0); + ncol=50; + found=kTRUE; + } + + if (paletteName.Contains("deepblue")) + { + Double_t s[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; + Double_t r[5] = { 0.00, 0.09, 0.18, 0.09, 0.00 }; + Double_t g[5] = { 0.01, 0.02, 0.39, 0.68, 0.97 }; + Double_t b[5] = { 0.17, 0.39, 0.62, 0.79, 0.97 }; + gStyle->CreateGradientColorTable(5, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("lightblue")) + { + Double_t s[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; + Double_t r[5] = { 0.00, 0.09, 0.18, 0.09, 0.00 }; + Double_t g[5] = { 0.00, 0.02, 0.40, 0.70, 1.00 }; + Double_t b[5] = { 0.00, 0.27, 0.51, 0.81, 1.00 }; + gStyle->CreateGradientColorTable(5, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("greyscale")) + { + double s[2] = {0.00, 1.00}; + double r[2] = {0.00, 1.00}; + double g[2] = {0.00, 1.00}; + double b[2] = {0.00, 1.00}; + gStyle->CreateGradientColorTable(2, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("glow1")) + { + double s[5] = {0., 0.10, 0.45, 0.75, 1.00}; + double r[5] = {0., 0.35, 0.85, 1.00, 1.00}; + double g[5] = {0., 0.10, 0.20, 0.73, 1.00}; + double b[5] = {0., 0.03, 0.06, 0.00, 1.00}; + gStyle->CreateGradientColorTable(5, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("glow2")) + { + double s[4] = {0.00, 0.50, 0.75, 1.00}; + double r[4] = {0.24, 0.67, 1.00, 1.00}; + double g[4] = {0.03, 0.04, 0.80, 1.00}; + double b[4] = {0.03, 0.04, 0.00, 1.00}; + gStyle->CreateGradientColorTable(4, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("glowsym")) + { + double s[8] = {0.00, 0.17, 0.39, 0.50, 0.55, 0.72, 0.88, 1.00}; + double r[8] = {0.09, 0.18, 0.09, 0.00, 0.35, 0.85, 1.00, 1.00}; + double g[8] = {0.70, 0.40, 0.02, 0.00, 0.10, 0.20, 0.73, 1.00}; + double b[8] = {0.81, 0.51, 0.27, 0.00, 0.03, 0.06, 0.00, 1.00}; + gStyle->CreateGradientColorTable(8, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("redish")) + { + double s[3] = {0., 0.5, 1.}; + double r[3] = {0., 1.0, 1.}; + double g[3] = {0., 0.0, 1.}; + double b[3] = {0., 0.0, 1.}; + gStyle->CreateGradientColorTable(3, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("bluish")) + { + double s[3] = {0., 0.5, 1.}; + double r[3] = {0., 0.0, 1.}; + double g[3] = {0., 0.0, 1.}; + double b[3] = {0., 1.0, 1.}; + gStyle->CreateGradientColorTable(3, s, r, g, b, ncol); + found=kTRUE; + } + + if (paletteName.Contains("small1")) + { + double s[4] = {0.00, 0.50, 0.95, 1.}; + double r[4] = {0.04, 0.28, 0.98, 1.}; + double g[4] = {0.28, 0.93, 0.03, 1.}; + double b[4] = {0.79, 0.11, 0.03, 1.}; + gStyle->CreateGradientColorTable(4, s, r, g, b, ncol); + found=kTRUE; + } + + if (inverse) + { + TArrayI c(ncol); + for (int i=0; iGetColorPalette(i); + gStyle->SetPalette(ncol, c.GetArray()); + } + + if (!found) + gLog << warn << "MH::SetPalette: Palette " << paletteName << " unknown... ignored." << endl; +} Index: trunk/MagicSoft/Mars/mhbase/MH.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.h (revision 7172) +++ trunk/MagicSoft/Mars/mhbase/MH.h (revision 7173) @@ -117,4 +117,6 @@ static TObject *FindObjectInPad(const char *name, TVirtualPad *pad=NULL); + static void SetPalette(TString paletteName, Int_t ncol=50); + ClassDef(MH, 2) //A base class for Mars histograms }; Index: trunk/MagicSoft/Mars/mhflux/MHDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7172) +++ trunk/MagicSoft/Mars/mhflux/MHDisp.cc (revision 7173) @@ -64,5 +64,5 @@ // MHDisp::MHDisp(const char *name, const char *title) - : fHilExt(0), fDisp(0) + : fHilExt(0), fDisp(0)//, fIsWobble(kFALSE) { // @@ -78,5 +78,5 @@ fHist.SetXTitle("x [\\circ]"); fHist.SetYTitle("y [\\circ]"); - fHist.SetZTitle("\\vartheta [deg^{2}]"); + fHist.SetZTitle("Eq.cts"); } @@ -130,12 +130,17 @@ return kFALSE; } + + //const Double_t xmax = fHist.GetXaxis()->GetXmax(); // Initialize all bins with a small (=0) value otherwise // these bins are not displayed for (int x=0; xGetBinCenter(x+1), - fHist.GetYaxis()->GetBinCenter(y+1), - 0.0, 1e-10); + for (int y=0; yGetBinCenter(x+1); + const Double_t cy = fHist.GetYaxis()->GetBinCenter(y+1); + //if (TMath::Hypot(cx, cy)RotationAngle(*fObservatory, *fTime); + // FIXME: Do wobble-rotation when subtracting? + if (!fHistOff/* && fIsWobble*/) + rho += TMath::Pi(); + // Get Disp from Parlist const Double_t disp = fDisp->GetVal(); - // Calculate where disp is pointing + // Calculate both postiions where disp could point TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta()); TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta()); - pos1 *= disp; - pos2 *= -disp; - + pos1 *= disp; // Vector of length disp in direction of shower + pos2 *= -disp; // Vector of length -disp in direction of shower + + // Move origin of vector to center-of-gravity of shower pos1 += hil->GetMean()*fMm2Deg; pos2 += hil->GetMean()*fMm2Deg; @@ -207,4 +217,5 @@ srcp = fSrcPos->GetXY(); + // Derotate all position around camera center by -rho if (rho!=0) { @@ -214,23 +225,34 @@ } + // Shift the source position to 0/0 pos1 -= srcp*fMm2Deg; pos2 -= srcp*fMm2Deg; - fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight); - fHist.Fill(pos2.X(), pos2.Y(), 0.0, w*(1-gweight)); + //const Double_t xmax = fHist.GetXaxis()->GetXmax(); + + // Fill histograms + //if (pos1.Mod()0.35) + sum += h.GetBinContent(x+1,y+1); + } + + return sum; +} void MHDisp::Paint(Option_t *o) @@ -240,7 +262,43 @@ pad->cd(1); + // Project on data onto yx-plane fHist.GetZaxis()->SetRange(0,9999); TH1 *h1=fHist.Project3D("yx_on"); - gStyle->SetPalette(1, 0); + + // Set Glow-palette (PaintSurface doesn't allow more than 99 colors) + MH::SetPalette(fHistOff?"glowsym":"glow1", 99); + h1->SetContour(99); + + Double_t scale = 1; + if (fHistOff) + { + // Project off data onto yx-plane and subtract it from on-data + fHistOff->GetZaxis()->SetRange(0,9999); + TH1 *h=fHistOff->Project3D("yx_off"); + + scale = -1; + + //if (!fIsWobble) + { + const Double_t h1off = GetOffSignal(*h1); + const Double_t hoff = GetOffSignal(*h); + scale = hoff==0?-1:-h1off/hoff; + } + + h1->Add(h, scale); + delete h; + + // Force calculation of minimum, maximum + h1->SetMinimum(); + h1->SetMaximum(); + + // Find maximum + const Double_t max = TMath::Max(TMath::Abs(h1->GetMinimum()), + TMath::Abs(h1->GetMaximum())); + + // Set new minimum, maximum + h1->SetMinimum(-max); + h1->SetMaximum( max); + } Int_t ix, iy, iz; @@ -304,5 +362,5 @@ const Double_t s = MMath::SignificanceLiMa(e, b); - h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f", x0, y0, func.GetParameter(2), s, e-b, b)); + h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f f=%.2f", x0, y0, func.GetParameter(2), s, e-b, b, scale)); } /* @@ -386,2 +444,38 @@ h->Draw(); } + +// -------------------------------------------------------------------------- +// +// The following resources are available: +// +// MHDisp.Wobble: on/off +// +/* +Int_t MHFalseSource::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "DistMin", print)) + { + rc = kTRUE; + SetMinDist(GetEnvValue(env, prefix, "DistMin", fMinDist)); + } + if (IsEnvDefined(env, prefix, "DistMax", print)) + { + rc = kTRUE; + SetMaxDist(GetEnvValue(env, prefix, "DistMax", fMaxDist)); + } + + if (IsEnvDefined(env, prefix, "DWMin", print)) + { + rc = kTRUE; + SetMinDW(GetEnvValue(env, prefix, "DWMin", fMinDist)); + } + if (IsEnvDefined(env, prefix, "DWMax", print)) + { + rc = kTRUE; + SetMaxDW(GetEnvValue(env, prefix, "DWMax", fMaxDist)); + } + + return rc; +} +*/ Index: trunk/MagicSoft/Mars/mhflux/MHDisp.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHDisp.h (revision 7172) +++ trunk/MagicSoft/Mars/mhflux/MHDisp.h (revision 7173) @@ -13,10 +13,15 @@ { private: - MHillasExt *fHilExt; //! - MParameterD *fDisp; //! + MHillasExt *fHilExt; //! + MParameterD *fDisp; //! - Double_t fM3lCut; //! - Double_t fXi; //! - Double_t fXiTheta; //! + Double_t fM3lCut; //! + Double_t fXi; //! + Double_t fXiTheta; //! + + //Bool_t fIsWobble; + + // MHDisp + Double_t GetOffSignal(TH1 &h) const; public: @@ -27,6 +32,10 @@ Bool_t Fill(const MParContainer *par, const Stat_t w=1); + // TObject void Paint(Option_t *o=""); void Draw(Option_t *o=""); + + // MParContainer + //Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); ClassDef(MHDisp, 1) //3D-histogram in alpha, x and y