source: trunk/MagicSoft/Mars/mhflux/MHThetaSqN.cc@ 7749

/* ======================================================================== *\
!
! *
! * 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 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2006
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MHThetaSqN
//
// For more detailes see MHAlpha.
//
// How to use this class:
//  - in your ganymed.rc add:
//      MJCut.NameHist: MHThetaSqN
//  - setup the number of off-source regions in your ganymed.rc by:
//      MHThetaSqN.NumOffSourcePos: 2
//  - switch on/off whether an off-theta cut should be done:
//      MHThetaSqN.DoOffCut: Yes,No
//  - and if necessary switch off the Theta cut in your Magic cuts:
//      Cut1.ThetaCut: None
//
//////////////////////////////////////////////////////////////////////////////
#include "MHThetaSqN.h"

#include <TVector2.h>

#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), fSignificanceCutLevel(1.7),
    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);

    fThetaSqCut = fSignificanceCutLevel*fFit.GetSignalIntegralMax()/1.7;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Return the value from fParemeter which should be filled into the plots
//
Double_t MHThetaSqN::GetVal() const
{
    return static_cast<const MParameterD*>(fParameter)->GetVal();
}

// --------------------------------------------------------------------------
//
// Abbreviation to set the value used by MHAlpha to fill the histogram
//
void MHThetaSqN::SetVal(Double_t val)
{
    static_cast<const MParameterD*>(fParameter)->SetVal(val);
}

// --------------------------------------------------------------------------
//
// Abbreviation to set the value used by MHAlpha to fill the histogram
//
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;

    // Calculate distance (theta-sq) to all (off-)source regions
    TArrayD dist(n);
    for (UInt_t i=0; i<n; i++)
    {
        const TVector2 src = const_cast<TVector2&>(src0).Rotate(i*rad);
        dist[i] = (src-org).Mod2();
    }

    // Processing off-data
    // Check if event's origin is in the on-regions
    if (!fOffData && fDoOffCut && dist[0]<fThetaSqCut)
        return kTRUE;

    for (UInt_t i=0; i<n; i++)
    {
        //   off: is in src region   on: is in off regions
        /// if (!fOffData && i==0) || (fOffData && i!=0)
        if ((bool)fOffData ^ (i==0) )
            continue;

        Stat_t w = weight;

        // Processing on-data
        if (fOffData && fDoOffCut)
        {
            /*
             static int cnt=0;
             if (dist[1+(cnt++%fNumOffSourcePos)]<fFit.GetSignalIntegralMax())
                continue;
             */

            // Check if event's origin is in one of the off-regions
            for (UInt_t j=1; j<n; j++)
                if (dist[j]<fThetaSqCut)
                {
                    w *= (float)(fNumOffSourcePos-1)/fNumOffSourcePos;
                    break;
                }
        }

        SetVal(dist[i]);

        if (!MHAlpha::Fill(NULL, w))
            return kFALSE;
    }

    if (!fOffData)
    {
        // off    0.4deg    0.3deg
        //  5: 0.200  ~86%   0.150
        //  4: 0.235  ~92%   0.176
        //  3: 0.283  ~96%   0.212
        //  2: 0.346  ~98%   0.260

        const Double_t dist = src0.Mod()*TMath::Sin(rad/2);
        if (dist<TMath::Sqrt(fFit.GetSignalIntegralMax()))
        {
            *fLog << warn << "WARNING - Source regions overlap: ";
            *fLog << "distance " << dist << " less than theta-sq cut ";
            *fLog << TMath::Sqrt(fFit.GetSignalIntegralMax()) << "!" << endl;
        }
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// 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.
//
// It takes fSkipHist* into account!
//
void MHThetaSqN::InitMapping(MHMatrix *mat, Int_t type)
{
    if (fMatrix)
        return;

    fMatrix = mat;

    fMap[0] = -1;
    fMap[1] = -1;
    fMap[2] = -1;
    fMap[3] = -1;
    fMap[4] = -1;

    if (!fSkipHistEnergy)
        if (type==0)
        {
            fMap[1] = fMatrix->AddColumn("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;
    }
    if (IsEnvDefined(env, prefix, "SignificanceCutLevel", print))
    {
        SetSignificanceCutLevel(GetEnvValue(env, prefix, "SignificanceCutLevel", fSignificanceCutLevel));
        rc = kTRUE;
    }
    return rc;
}