Index: trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc (revision 7716) +++ trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc (revision 7716) @@ -0,0 +1,326 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 5/2006 +! +! Copyright: MAGIC Software Development, 2000-2006 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHThetaSqN +// +// For more detailes see MHAlpha. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHThetaSqN.h" + +#include + +#include "MParameters.h" +#include "MSrcPosCam.h" +#include "MGeomCam.h" +#include "MHMatrix.h" +#include "MHillas.h" + +#include "MBinning.h" +#include "MParList.h" +#include "MTaskList.h" +#include "MParameters.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHThetaSqN); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor +// +MHThetaSqN::MHThetaSqN(const char *name, const char *title) + : MHAlpha(name, title), fDisp(0), fSrcPosCam(0), + fNumBinsSignal(3), fNumBinsTotal(75), fNumOffSourcePos(3), fDoOffCut(kTRUE) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHThetaSqN"; + fTitle = title ? title : "Theta Squared plot"; + + fNameParameter = "ThetaSquared"; + + fHist.SetName("Theta"); + fHist.SetTitle("Theta"); + fHist.SetZTitle("\\vartheta^{2} [deg^{2}]"); + fHist.SetDirectory(NULL); + + // Main histogram + fHistTime.SetName("Theta"); + fHistTime.SetXTitle("\\vartheta^{2} [deg^{2}]"); + fHistTime.SetDirectory(NULL); + + MBinning binsa, binse, binst; + binsa.SetEdges(75, 0, 1.5); + //binsa.SetEdges(arr); + binse.SetEdgesLog(15, 10, 100000); + binst.SetEdgesASin(67, -0.005, 0.665); + binsa.Apply(fHistTime); + + MH::SetBinning(&fHist, &binst, &binse, &binsa); + + fParameter = new MParameterD; +} + +MHThetaSqN::~MHThetaSqN() +{ + delete fParameter; +} + +// -------------------------------------------------------------------------- +// +// Overwrites the binning in Alpha (ThetaSq) with a binning for which +// the upper edge of the 5th bin (bin=5) fit the signal integration window. +// In total 75 bins are setup. +// +// In case of fOffData!=NULL the binnings are taken later from fOffData anyhow. +// +Bool_t MHThetaSqN::SetupFill(const MParList *pl) +{ + // Default is from default fitter + // if a user defined fitter is in the parlist: use this range + MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter"); + if (!fit) + fit = &fFit; + + MParameterD *cut = (MParameterD*)pl->FindObject("ThetaSquaredCut", "MParameterD"); + if (cut) + fit->SetSignalIntegralMax(cut->GetVal()); + + // Get Histogram binnings + MBinning binst, binse; + binst.SetEdges(fHist, 'x'); + binse.SetEdges(fHist, 'y'); + + // Calculate bining which fits alpha-cut + const Double_t intmax = fit->GetSignalIntegralMax(); + const UInt_t nbins = fNumBinsTotal; + const UInt_t nsig = fNumBinsSignal; + + MBinning binsa(nbins, 0, nbins*intmax/nsig); + + // Apply binning + binsa.Apply(fHistTime); + MH::SetBinning(&fHist, &binst, &binse, &binsa); + + // Remark: Binnings might be overwritten in MHAlpha::SetupFill + if (!MHAlpha::SetupFill(pl)) + return kFALSE; + + fDisp = (MParameterD*)pl->FindObject("Disp", "MParameterD"); + if (!fDisp) + { + *fLog << err << "Disp [MParameterD] not found... abort." << endl; + return kFALSE; + } + fHillas = (MHillas*)pl->FindObject("MHillas"); + if (!fHillas) + { + *fLog << err << "MHillas not found... abort." << endl; + return kFALSE; + } + fSrcPosCam = (MSrcPosCam*)pl->FindObject("MSrcPosCam"); + if (!fSrcPosCam) + { + *fLog << err << "MSrcPosCam not found... abort." << endl; + return kFALSE; + } + + MGeomCam *geom = (MGeomCam*)pl->FindObject("MGeomCam"); + if (!geom) + { + *fLog << err << "MGeomCam not found... abort." << endl; + return kFALSE; + } + fMm2Deg = geom->GetConvMm2Deg(); + + if (fFit.GetScaleMode()==MAlphaFitter::kNone) + fFit.SetScaleUser(1./fNumOffSourcePos); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Return the value from fParemeter which should be filled into the plots +// +Double_t MHThetaSqN::GetVal() const +{ + return static_cast(fParameter)->GetVal(); +} + +void MHThetaSqN::SetVal(Double_t val) +{ + static_cast(fParameter)->SetVal(val); +} + +TVector2 MHThetaSqN::GetVec(const TVector2 &v, Int_t n1) const +{ + if (!fMatrix) + return v; + + return TVector2( (*fMatrix)[n1], (*fMatrix)[n1+1] ); +} + +Bool_t MHThetaSqN::Fill(const MParContainer *par, const Stat_t weight) +{ + const TVector2 mean(GetVec(fHillas->GetMean(), 6)); + const TVector2 norm(GetVec(fHillas->GetNormAxis(), 8)); + + const TVector2 org = mean*fMm2Deg + norm*fDisp->GetVal(); + + // In the case we are filling off-data src0 contains the anti-source position + TVector2 src0(GetVec(fSrcPosCam->GetXY(), 10)*fMm2Deg); + if (!fOffData) + src0 *= -1; + + const UInt_t n = fNumOffSourcePos+1; + const Float_t rad = TMath::TwoPi()/n; + + for (UInt_t i=0; i(src0).Rotate(i*rad); + const Double_t d = (src-org).Mod2(); + + // off: is in src region on: is in off regions + if ((bool)fOffData ^ (i==0) /*(!fOffData && i==0) || (fOffData && i!=0)*/) + { + if (dAddColumn("MEnergyEst.fVal"); + fMap[2] = -1; + } + else + { + fMap[1] = -1; + fMap[2] = fMatrix->AddColumn("MHillas.fSize"); + } + + if (!fSkipHistTheta) + fMap[3] = fMatrix->AddColumn("MPointingPos.fZd"); + + fMap[5] = fMatrix->AddColumn("Disp.fVal"); + fMap[6] = fMatrix->AddColumn("MHillas.fCosDelta"); + fMap[7] = fMatrix->AddColumn("MHillas.fSinDelta"); + fMap[8] = fMatrix->AddColumn("MHillas.fMeanX"); + fMap[9] = fMatrix->AddColumn("MHillas.fMeanY"); + fMap[10] = fMatrix->AddColumn("MSrcPosCam.fX"); + fMap[11] = fMatrix->AddColumn("MSrcPosCan.fY"); + + // if (!fSkipHistTime) + // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime"); +} + +Int_t MHThetaSqN::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Int_t rc = MHAlpha::ReadEnv(env, prefix, print); + if (rc==kERROR) + return kERROR; + + if (IsEnvDefined(env, prefix, "NumBinsSignal", print)) + { + SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal)); + rc = kTRUE; + } + if (IsEnvDefined(env, prefix, "NumBinsTotal", print)) + { + SetNumBinsTotal(GetEnvValue(env, prefix, "NumBinsTotal", (Int_t)fNumBinsTotal)); + rc = kTRUE; + } + if (IsEnvDefined(env, prefix, "NumOffSourcePos", print)) + { + SetNumOffSourcePos(GetEnvValue(env, prefix, "NumOffSourcePos", (Int_t)fNumOffSourcePos)); + rc = kTRUE; + } + if (IsEnvDefined(env, prefix, "DoOffCut", print)) + { + SetDoOffCut(GetEnvValue(env, prefix, "DoOffCut", fDoOffCut)); + rc = kTRUE; + } + return rc; +}