Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 2224) +++ /trunk/MagicSoft/Mars/Changelog (revision 2225) @@ -1,3 +1,18 @@ -*-*- END OF LINE -*-*- + + 2003/06/24: Thomas Bretz + + * manalysis/MCT1SupercutsCalc.[h,cc]: + - implemented Mapping for Supercuts + - changed data member arrays to TArrayD + + * manalysis/MEnergyEstParam.h: + - added a comment + + * mhist/MHHadronness.[h,cc]: + - implemented mapping + - implemented calculating Acc_g/sqrt(Acc_h) for filtercuts + + 2003/06/23: Thomas Bretz Index: /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2224) +++ /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2225) @@ -33,7 +33,4 @@ #include -#include "MLog.h" -#include "MLogManip.h" - #include "MParList.h" #include "MHillasExt.h" @@ -43,4 +40,8 @@ #include "MGeomCam.h" #include "MHadronness.h" +#include "MHMatrix.h" + +#include "MLog.h" +#include "MLogManip.h" ClassImp(MCT1SupercutsCalc); @@ -50,4 +51,14 @@ void MCT1SupercutsCalc::InitParams() { + fLengthUp.Set(8); + fLengthLo.Set(8); + fWidthUp.Set(8); + fWidthLo.Set(8); + fDistUp.Set(8); + fDistLo.Set(8); + fAsymUp.Set(8); + fAsymLo.Set(8); + fAlphaUp.Set(8); + //--------------------------------- // cut parameters @@ -151,4 +162,23 @@ Int_t MCT1SupercutsCalc::PreProcess(MParList *pList) { + MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!cam) + { + *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; + return kFALSE; + } + + fMm2Deg = cam->GetConvMm2Deg(); + + fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName); + if (!fHadronness) + { + *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl; + return kFALSE; + } + + if (fMatrix) + return kTRUE; + fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); if (!fHil) @@ -172,21 +202,4 @@ } - MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!cam) - { - *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; - return kFALSE; - } - - fMm2Deg = cam->GetConvMm2Deg(); - - fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName); - if (!fHadronness) - { - *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl; - return kFALSE; - } - - return kTRUE; } @@ -196,5 +209,9 @@ // Calculation of upper and lower limits // -Double_t MCT1SupercutsCalc::CtsMCut(Double_t *a, Double_t ls, Double_t ct, +Double_t MCT1SupercutsCalc::CtsMCut( +#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00) +const +#endif + TArrayD &a, Double_t ls, Double_t ct, Double_t ls2, Double_t dd2) { @@ -209,5 +226,4 @@ // ct: Cos(ZA.) - dNOMCOSZA // dd2: DIST^2 - const Double_t limit = a[0] + a[1] * dd2 + a[2] * ct + @@ -225,4 +241,38 @@ return limit; +} + +// -------------------------------------------------------------------------- +// +// Returns the mapped value from the Matrix +// +Double_t MCT1SupercutsCalc::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 MCT1SupercutsCalc::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + fMatrix = mat; + + fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta"); + fMap[1] = fMatrix->AddColumn("MHillas.fWidth"); + fMap[2] = fMatrix->AddColumn("MHillas.fLength"); + fMap[3] = fMatrix->AddColumn("MHillas.fSize"); + fMap[4] = fMatrix->AddColumn("MHillas.fMeanX"); + fMap[5] = fMatrix->AddColumn("MHillas.fMeanY"); + fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist"); } @@ -240,26 +290,25 @@ const Double_t kNomCosZA = 1.0; - const Double_t newdist = fHilSrc->GetDist() * fMm2Deg; - - const Double_t xbar = fHil->GetMeanX(); - const Double_t ybar = fHil->GetMeanY(); - const Double_t dist = sqrt(xbar*xbar + ybar*ybar) * fMm2Deg; - const Double_t dd2 = dist * dist; - - const Double_t dmls = log(fHil->GetSize()) - kNomLogSize; + const Double_t theta = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta(); + const Double_t width0 = fMatrix ? GetVal(1) : fHil->GetWidth(); + const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength(); + const Double_t size = fMatrix ? GetVal(3) : fHil->GetSize(); + const Double_t meanx = fMatrix ? GetVal(4) : fHil->GetMeanX(); + const Double_t meany = fMatrix ? GetVal(5) : fHil->GetMeanY(); + const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist(); + + const Double_t newdist = dist0 * fMm2Deg; + + const Double_t dist2 = meanx*meanx + meany*meany; + const Double_t dd2 = dist2*fMm2Deg; + const Double_t dist = sqrt(dist2) * fMm2Deg; + + const Double_t dmls = log(size) - kNomLogSize; const Double_t dmls2 = dmls * dmls; - const Double_t dmcza = cos(fMcEvt->GetTelescopeTheta()) - kNomCosZA; - - const Double_t length = fHil->GetLength() * fMm2Deg; - const Double_t width = fHil->GetWidth() * fMm2Deg; - - //*fLog << "MCT1SupercutsCalc::Process; dmls, dmcza, dmls2, dd2 = " - // << dmls << ", " << dmcza << ", " << dmls2 << ", " - // << dd2 << endl; - - //*fLog << "MCT1SupercutsCalc::Process; newdist, dist, length, width = " - // << newdist << ", " << dist << ", " << length << ", " - // << width << endl; + const Double_t dmcza = cos(theta) - kNomCosZA; + + const Double_t length = length0 * fMm2Deg; + const Double_t width = width0 * fMm2Deg; if (newdist < 1.05 && Index: /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2224) +++ /trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2225) @@ -2,12 +2,10 @@ #define MARS_MCT1SupercutsCalc -///////////////////////////////////////////////////////////////////////////// -// // -// MCT1SupercutsCalc // -// // -///////////////////////////////////////////////////////////////////////////// - #ifndef MARS_MFilter #include "MFilter.h" +#endif + +#ifndef ROOT_TArrayD +#include #endif @@ -19,4 +17,5 @@ class MGeomCam; class MHadronness; +class MHMatrix; @@ -35,15 +34,18 @@ Double_t fMm2Deg; + Int_t fMap[7]; + MHMatrix *fMatrix; + //--------------------------------- // cut parameters - Double_t fLengthUp[8]; - Double_t fWidthUp[8]; - Double_t fDistUp[8]; - Double_t fLengthLo[8]; - Double_t fWidthLo[8]; - Double_t fDistLo[8]; - Double_t fAsymUp[8]; - Double_t fAsymLo[8]; - Double_t fAlphaUp[8]; + TArrayD fLengthUp; + TArrayD fLengthLo; + TArrayD fWidthUp; + TArrayD fWidthLo; + TArrayD fDistUp; + TArrayD fDistLo; + TArrayD fAsymUp; + TArrayD fAsymLo; + TArrayD fAlphaUp; //--------------------------------- @@ -53,4 +55,14 @@ Int_t Process(); + void Set(TArrayD &a, const TArrayD &b) { if (a.GetSize()==b.GetSize()) a=b; } + Double_t GetVal(Int_t i) const; + + Double_t CtsMCut( +#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00) +const +#endif + TArrayD &a, Double_t ls, Double_t ct, + Double_t ls2, Double_t dd2); + public: MCT1SupercutsCalc(const char *hilname="MHillas", @@ -58,9 +70,29 @@ const char *name=NULL, const char *title=NULL); - Double_t CtsMCut(Double_t *a, Double_t ls, Double_t ct, - Double_t ls2, Double_t dd2); - void SetHadronnessName(const TString name) { fHadronnessName = name; } TString GetHadronnessName() const { return fHadronnessName; } + + void InitMapping(MHMatrix *mat); + void StopMapping() { InitMapping(NULL); } + + void SetLengthUp(const TArrayD &d) { Set(fLengthUp, d); } + void SetLengthLo(const TArrayD &d) { Set(fLengthLo, d); } + void SetWidthUp (const TArrayD &d) { Set(fWidthUp, d); } + void SetWidthLo (const TArrayD &d) { Set(fWidthLo, d); } + void SetDistUp (const TArrayD &d) { Set(fDistUp, d); } + void SetDistLo (const TArrayD &d) { Set(fDistLo, d); } + void SetAsymUp (const TArrayD &d) { Set(fAsymUp, d); } + void SetAsymLo (const TArrayD &d) { Set(fAsymLo, d); } + void SetAlphaUp (const TArrayD &d) { Set(fAlphaUp, d); } + + const TArrayD &GetLengthUp() const { return fLengthUp; } + const TArrayD &GetLengthLo() const { return fLengthLo; } + const TArrayD &GetWidthUp() const { return fWidthUp; } + const TArrayD &GetWithLo() const { return fWidthLo; } + const TArrayD &GetDistUp() const { return fDistUp; } + const TArrayD &GetDistLo() const { return fDistLo; } + const TArrayD &GetAsymUp() const { return fAsymUp; } + const TArrayD &GetAsymLo() const { return fAsymLo; } + const TArrayD &GetAlphaUp() const { return fAlphaUp; } ClassDef(MCT1SupercutsCalc, 0) // A class to evaluate the Supercuts Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h (revision 2224) +++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h (revision 2225) @@ -19,5 +19,5 @@ { private: - Int_t fMap[100]; + Int_t fMap[100]; // FIXME! MHMatrix *fMatrix; Index: /trunk/MagicSoft/Mars/mhist/MHHadronness.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 2224) +++ /trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 2225) @@ -56,4 +56,9 @@ // MFillH. // +// If you are using filtercuts which gives you only two discrete values +// of the hadronness (0.25 and 0.75) you may want to use MHHadronness(2). +// If you request Q05() from such a MHHadronness instance you will get +// Acc_g/sqrt(Acc_h) +// //////////////////////////////////////////////////////////////////////////// #include "MHHadronness.h" @@ -65,11 +70,12 @@ #include +#include "MParList.h" +#include "MBinning.h" +#include "MHMatrix.h" +#include "MHadronness.h" + #include "MLog.h" #include "MLogManip.h" -#include "MParList.h" -#include "MBinning.h" -#include "MHadronness.h" - #include "MMcEvt.hxx" @@ -84,4 +90,5 @@ // MHHadronness::MHHadronness(Int_t nbins, const char *name, const char *title) + : fMatrix(NULL) { // @@ -95,16 +102,16 @@ fGraph->SetMarkerStyle(kFullDotSmall); - fGhness = new TH1D("Ghness", "H. Gammas", nbins, 0, 1); - fPhness = new TH1D("Phness", "H. Hadrons", nbins, 0, 1); + fGhness = new TH1D("Ghness", "Acceptance vs. Hadronness (Gammas)", nbins, 0, 1); + fPhness = new TH1D("Phness", "Acceptance vs. Hadronness (Hadrons)", nbins, 0, 1); fGhness->SetXTitle("Hadronness"); fPhness->SetXTitle("Hadronness"); - fGhness->SetYTitle("Counts"); - fPhness->SetYTitle("Counts"); + fGhness->SetYTitle("Acceptance"); + fPhness->SetYTitle("Acceptance"); fPhness->SetLineColor(kRed); fGhness->SetDirectory(NULL); fPhness->SetDirectory(NULL); - fIntGhness = new TH1D("AccGammas", "Acceptance", nbins, 0, 1); - fIntPhness = new TH1D("AccHadrons", "Acceptance", nbins, 0, 1); + fIntGhness = new TH1D("AccGammas", "Integral Acceptance vs. Hadronness (Gammas)", nbins, 0, 1); + fIntPhness = new TH1D("AccHadrons", "Integral Acceptance vs. Hadronness (Hadrons)", nbins, 0, 1); fIntGhness->SetXTitle("Hadronness"); fIntPhness->SetXTitle("Hadronness"); @@ -145,12 +152,18 @@ Bool_t MHHadronness::SetupFill(const MParList *plist) { - fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; - return kFALSE; + if (!fMatrix) + { + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + return kFALSE; + } } fHadronness = (MHadronness*)plist->FindObject("MHadronness"); + + fGhness->Reset(); + fPhness->Reset(); /* @@ -200,5 +213,7 @@ return kCONTINUE; - if (fMcEvt->GetPartId()==kGAMMA) + const Int_t particleid = fMatrix ? (Int_t)(*fMatrix)[fMap] : fMcEvt->GetPartId(); + + if (particleid==kGAMMA) fGhness->Fill(h, w); else @@ -208,7 +223,40 @@ } +// -------------------------------------------------------------------------- +// +// Returns the quality factor at gamma acceptance 0.5. +// +// If the histogram containes only two bins we return the +// naive quality: ag/sqrt(ah) +// with ag the acceptance of gammas and ah the acceptance of hadrons. +// +// You can use this (nbins=2) in case of some kind of filter cuts giving +// only a result: gamma yes/no (means discrete values of hadronness 0.25 +// or 0.75) +// +// FIXME: In the later case weights cannot be used! +// Float_t MHHadronness::GetQ05() const { - Int_t n = fQfac->GetN(); + if (fGhness->GetNbinsX()==2) + { + // acceptance of all gamma-like gammas (h<0.5) + const Double_t ig = fGhness->GetBinContent(1); + + // acceptance of all gamma-like hadrons (h<0.5) + const Double_t ip = fPhness->GetBinContent(1); + + if (ip==0) + return 0; // FIXME! + + // naive quality factor + const Double_t q = ig / sqrt(ip); + + *fLog << all << ip << "/" << ig << ": " << q << endl; + + return q; + } + + const Int_t n = fQfac->GetN(); Double_t val1x=0; @@ -479,2 +527,26 @@ } } + +// -------------------------------------------------------------------------- +// +// You can use this function if you want to use a MHMatrix instead of +// MMcEvt. 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) MHHadronness::Fill +// will take the values from the matrix instead of the containers. +// +void MHHadronness::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + fMatrix = mat; + fMap = fMatrix->AddColumn("MMcEvt.fPartId"); +} + +void MHHadronness::StopMapping() +{ + fMatrix = NULL; +} + Index: /trunk/MagicSoft/Mars/mhist/MHHadronness.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHadronness.h (revision 2224) +++ /trunk/MagicSoft/Mars/mhist/MHHadronness.h (revision 2225) @@ -11,4 +11,5 @@ class MMcEvt; class MHadronness; +class MHMatrix; class MHHadronness : public MH @@ -17,4 +18,6 @@ const MMcEvt *fMcEvt; //! const MHadronness *fHadronness; //! + MHMatrix *fMatrix; //! + Int_t fMap; //! TH1D* fPhness; //-> Hadrons Hadronness @@ -48,4 +51,7 @@ Bool_t Finalize(); + void InitMapping(MHMatrix *mat); + void StopMapping(); + void Print(Option_t *option="") const;