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

/* ======================================================================== *\
!
! *
! * 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 "MMcTrig.hxx" 
#include "MHMcRate.h"

ClassImp(MMcTriggerRateCalc)

// --------------------------------------------------------------------------
//
//  default constructor
//  set member data to initial values
//
MMcTriggerRateCalc::MMcTriggerRateCalc (const char *name, const char *title)
{
  //  default constructor
  *fName  = name  ? name  : "MMcTriggerRateCalc";
  *fTitle = title ? title : "Task to calc the trigger rate";

  fDimension=0;
  for (int i=0;i<10;i++){
  fTrigger[i]=0;
  }
  fShowers = 0;
  fAnalShow = 100000;

  fPartId=14;

  fSpecInd=2.75;
  fFlux0=10.92e-2;
} 

// --------------------------------------------------------------------------
//
//  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 (int dim, int part, float *trigbg,
                                float simbg, float rate,
                                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]=trigbg[i];
  }

  fShowers = 0;
  fAnalShow = simbg;

  fPartId=part;

  fSpecInd=rate;
  fFlux0=0;
}

// --------------------------------------------------------------------------
//
//  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: sSpecInd
//      flux0; fFlux0
//
MMcTriggerRateCalc::MMcTriggerRateCalc (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]=trigbg[i];
  }

  fShowers = 0;
  fAnalShow = simbg;

  fPartId=part;

  fSpecInd=spec;
  fFlux0=flux0;
} 

// --------------------------------------------------------------------------
//
//  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
  int i;               // bucle variable
  char auxname[15];    // string to write container names

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

  if(fDimension>0)
    for(i=0;i<fDimension;i++){
      sprintf(auxname,"MMcTrig;%i",i+1);
      fMcTrig[i] = (MMcTrig*)pList->FindObject(auxname);
      if (!fMcTrig[i])
        {
          *fLog << dbginf << "MMcTrig not found... aborting." << endl;
          return kFALSE;
        } 
      sprintf(auxname,"MHMcRate;%i",i+1);
      fMcRate[i] = (MHMcRate*)pList->FindCreateObj(auxname);
      if (!fMcTrig[i])
        {
          return kFALSE;
        } 
      fMcRate[i]->SetParticle(fPartId);
      fMcRate[i]->SetBackground(fTrigger[i],fAnalShow);
      fTrigger[i]=0;
    }
  else if (fDimension==0)
    {  
      fMcTrig[0] = (MMcTrig*)pList->FindObject("MMcTrig");
      if (!fMcTrig[0])
        {
          *fLog << dbginf << "MMcTrig not found... aborting." << endl;
          return kFALSE;
        } 
      fMcRate[0] = (MHMcRate*)pList->FindCreateObj("MHMcRate");
      if (!fMcTrig[0])
        {
          return kFALSE;
        } 
      fMcRate[0]->SetParticle(fPartId);
      fMcRate[0]->SetBackground(fTrigger[0],fAnalShow);
      fTrigger[0]=0;
      fDimension=1;
    }   
  else
     {
      sprintf(auxname,"MMcTrig;%i.",-fDimension);
      fMcTrig[0] = (MMcTrig*)pList->FindObject(auxname);
      if (!fMcTrig[0])
        {
          *fLog << dbginf << "MMcTrig not found... aborting." << endl;
          return kFALSE;
        } 
     sprintf(auxname,"MHMcRate;%i.",-fDimension);
     fMcRate[0] = (MHMcRate*)pList->FindCreateObj(auxname);
      if (!fMcTrig[0])
        {
          return kFALSE;
        } 
      fMcRate[0]->SetParticle(fPartId);
      fMcRate[0]->SetBackground(fTrigger[0],fAnalShow);
      fTrigger[0]=0;
      fDimension=1;
    }
  
  fAnalShow=0.0;
  *fLog<<dbginf<< "Pre Process done"<<endl;
  
  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 () 
{ 
  Float_t Ener, Theta, Phi, Param;

  //  Counting analysed and simulated showers

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

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

  Theta=fMcEvt->GetTheta();
  Phi=fMcEvt->GetPhi();
  Param=fMcEvt->GetImpact();
  Ener=fMcEvt->GetEnergy()/1000.0;

  //  Counting number of triggers

  for(int i=0;i<fDimension;i++){
  fTrigger[i]=fTrigger[i]+fMcTrig[i]->GetFirstLevel();
  fMcRate[i]->Boundaries(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++){
    fMcRate[i]->CalcRate(fTrigger[i],fAnalShow,fShowers);
    fMcRate[i]->Print();
    fMcRate[i]->Draw("S");
  }

  return kTRUE ; 
}