source: trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc@ 959

/* ======================================================================== *\
!
! *
! * 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): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!   Author(s): Thomas Bretz  12/2000 (tbretz@uni-sw.gwdg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

#include "MMcTriggerRateCalc.h"

#include "MLog.h"
#include "MLogManip.h"
#include "MParList.h"

#include "MHMcRate.h"
#include "MMcEvt.hxx"
#include "MMcTrig.hxx" 

ClassImp(MMcTriggerRateCalc);

void MMcTriggerRateCalc::Init(int dim, int part, float *trigbg,
                              float simbg, float spec, float flux0,
                              const char *name, const char *title)
{
    *fName  = name  ? name  : "MMcTriggerRateCalc";
    *fTitle = title ? title : "Task to calc the trigger rate ";

    fDimension=dim;

    for (int i=0;i<10;i++)
        fTrigger[i] = dim&&trigbg ? trigbg[i] : 0;

    fShowers = 0;
    fAnalShow = simbg;

    fPartId=part;

    fSpecInd=spec;
    fFlux0=flux0;
}

// --------------------------------------------------------------------------
//
//  overloaded constructor I
//
//      dim: fDimension
//      part: fPartId
//      *trigbg: number of shower from bacground that triggers
//               a given trigger condition.
//      simbg: Number of simulated showers for the bacground
//      rate: rate of incident showers
//

MMcTriggerRateCalc::MMcTriggerRateCalc(float rate, int dim, int part,
                                       float *trigbg, float simbg,
                                       const char *name, const char *title)
{
    Init(dim, part, trigbg, simbg, rate, 0, name, title);
}


// --------------------------------------------------------------------------
//
//  overloaded constructor II
//
//      dim: fDimension
//      part: fPartId
//      *trigbg: number of shower from bacground that triggers
//               a given trigger condition.
//      simbg: Number of simulated showers for the bacground
//      spec: spectral index
//      flux0; fFlux0
//

MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, int part, float *trigbg,
                                       float simbg, float spec, float flux0,
                                       const char *name, const char *title)
{
    Init(dim, part, trigbg, simbg, spec, flux0, name, title);
}


// --------------------------------------------------------------------------
//
//  The PreProcess connects the raw data with this task. It checks if the 
//  input containers exist, if not a kFalse flag is returned. It also checks
//  if the output contaniers exist, if not they are created.
//  This task can read either Montecarlo files with multiple trigger
//  options, either Montecarlo files with a single trigger option.
//
Bool_t MMcTriggerRateCalc::PreProcess (MParList *pList)
{
    // connect the raw data with this task

    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    if (!fMcEvt)
    {
        *fLog << dbginf << "MMcEvt not found... aborting." << endl;
        return kFALSE;
    }

    const UInt_t from = fDimension<=0 ? -fDimension : 0;

    if (fDimension<=0)
        fDimension = -fDimension;

    fMcTrig = new TObjArray(pList->FindObjectList("MMcTrig", from, fDimension));
    if (fMcTrig->GetEntriesFast() != fDimension)
    {
        *fLog << dbginf << "Error: Not all requested MMcTrig objects are available...aborting." << endl;
        return kFALSE;
    }

    fMcRate = new TObjArray(pList->FindObjectList("MHMcRate", from, fDimension));
    if (fMcRate->GetEntriesFast() != fDimension)
    {
        *fLog << dbginf << "Error: Not all requested MHMcRate objects are available...aborting." << endl;
        return kFALSE;
    }

    for (int i=0; i<fDimension; i++)
    {
        MHMcRate &rate = *GetRate(i);

        rate.SetParticle(fPartId);
        rate.SetBackground(fTrigger[i], fAnalShow);

        fTrigger[i]=0;
    }

    fAnalShow=0.0;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  The Process-function counts the number of simulated showers, the
//  number of analised showers and the number of triggers. It also updates
//  the limits for theta, phi, energy and impact parameter in the
//  MHMcRate container.
//
Bool_t MMcTriggerRateCalc::Process()
{
    //  Counting analysed and simulated showers

    fShowers++;
    if (fMcEvt->GetPhotElfromShower())
        fAnalShow++;

    //  Getting angles, energy and impact parameter to set boundaries

    const Float_t theta=fMcEvt->GetTheta();
    const Float_t phi  =fMcEvt->GetPhi();
    const Float_t param=fMcEvt->GetImpact();
    const Float_t ener =fMcEvt->GetEnergy()/1000.0;

    //  Counting number of triggers

    for (int i=0; i<fDimension; i++)
    {
        fTrigger[i] += GetTrig(i)->GetFirstLevel();

        GetRate(i)->UpdateBoundaries(ener, theta, phi, param);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  The PostProcess-function calculates and shows the trigger rate
//
Bool_t MMcTriggerRateCalc::PostProcess()
{
    // Computing trigger rate and showing it
    for (int i=0; i<fDimension; i++)
    {
        MHMcRate &rate = *GetRate(i);

        rate.CalcRate(fTrigger[i], fAnalShow, fShowers);
        // rate.Print();
        // rate.Draw("S");
    }

    return kTRUE;
}