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

/* ======================================================================== *\
!
! *
! * 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  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
!   Author(s): Harald Kornmayer 1/2001
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!   Modified 4/7/2002 Abelardo Moralejo: now the dimension of fTrigger is
!    set dinamically, to allow an arbitrary large number of trigger 
!    conditions to be processed.
!
!
\* ======================================================================== */

#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,
                              const char *name, const char *title)
{
    fName  = name  ? name  : "MMcTriggerRateCalc";
    fTitle = title ? title : "Task to calc the trigger rate ";

    fMcTrig = NULL;
    fMcRate = NULL;

    fShowers = 0;
    fAnalShow = simbg;

    fPartId=part;

    fFirst = dim>0 ?   1 : -dim;
    fLast  = dim>0 ? dim : -dim;

    fNum = fLast-fFirst+1;

    fTrigger = new float[fNum];

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

    AddToBranchList("MMcEvt.fImpact");
    AddToBranchList("MMcEvt.fEnergy");
    AddToBranchList("MMcEvt.fPhi");
    AddToBranchList("MMcEvt.fTheta");
    AddToBranchList("MMcEvt.fPhotElfromShower");
    AddToBranchList("MMcTrig", "fNumFirstLevel", fFirst, fLast);
}

// --------------------------------------------------------------------------
//
//  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 background
//      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, 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 background
//
MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, int part, float *trigbg,
                                       float simbg,
                                       const char *name, const char *title)
{
    Init(dim, part, trigbg, simbg, name, title);
}

MMcTriggerRateCalc::~MMcTriggerRateCalc()
{
    if (fMcTrig)
        delete fMcTrig;

    if (fMcRate)
        delete fMcRate;
}


// --------------------------------------------------------------------------
//
//  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 << err << dbginf << "MMcEvt not found... aborting." << endl;
        return kFALSE;
    }

    UInt_t num;

    fMcTrig = new TObjArray(pList->FindObjectList("MMcTrig", fFirst, fLast));
    num = fMcTrig->GetEntriesFast();
    if (num != fNum)
    {
        *fLog << err << dbginf << fNum << " MMcTrig objects requested, ";
        *fLog << num << " are available... aborting." << endl;
        return kFALSE;
    }

    fMcRate = new TObjArray(pList->FindObjectList("MHMcRate", fFirst, fLast));
    num = fMcRate->GetEntriesFast();
    if (num != fNum)
    {
        *fLog << err << dbginf << fNum << " MHMcRate objects requested, ";
        *fLog << num << " are available... aborting." << endl;
        return kFALSE;
    }

    for (UInt_t i=0; i<fNum; i++)
    {
        MHMcRate &rate = *GetRate(i);

        rate.SetParticle(fPartId);
        switch (fPartId)
        {
        case kPROTON:
            rate.SetFlux(0.1091, 2.75);
            break;
        case kHELIUM:
            rate.SetFlux(0.0660, 2.62);
            break;
        default:
            *fLog << err << dbginf << "Unknown incident flux parameters for ";
            *fLog << fPartId<< " particle Id ... aborting." << endl;
            return kFALSE;
        }
        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 (UInt_t i=0; i<fNum; i++)
    {
        if (GetTrig(i)->GetFirstLevel())
            fTrigger[i] ++; 

        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
    //
    for (UInt_t i=0; i<fNum; i++)
        GetRate(i)->CalcRate(fTrigger[i], fAnalShow, fShowers);

    return kTRUE;
}