source: trunk/MagicSoft/Mars/mhist/MHMcRate.cc@ 887

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

#include "MLog.h"

ClassImp(MHMcRate)

// --------------------------------------------------------------------------
//
//  default constructor
//  fills all member data with initial values
//
MHMcRate::MHMcRate (const char *name, const char *title)
{
  *fName  = name  ? name  : "MMcTriggerRate";
  *fTitle = title ? title : "Task to calc the collection area ";

  fPartId=0;              // Type of particle
  
  fEnergyMax=0.0;         // Maximum Energy in GeV
  fEnergyMin=1000000.0;   // Minimum Energy in GeV
  
  fThetaMax=0.0;          // Maximum theta angle of run 
  fThetaMin=370.0;        // Minimum theta angle of run
  fPhiMax=0.0;            // Maximum phi angle of run 
  fPhiMin=370.0;          // Minimum phi angle of run
  
  fImpactMax=0.0;         // Maximum impact parameter
  fImpactMin=100000.0;    // Minimum impact parameter
  
  fBackTrig=-1.0;         // Number of triggers from background
  fBackSim=-1.0;          // Number of simulated showers for the background
  
  fSpecIndex=0.0;         // dn/dE = k * e^{- fSpecIndex}
  fFlux0=-1.0;            // dn/dE = fFlux0 * E^{-a}
  
  fShowRate= -1.0;        // Showers rate in Hz
  fShowRateError=0.0;     // Estimated error of shower rate in Hz
  
  fTrigRate= -1.0;        // Trigger rate in Hz
  fTrigRateError= -1.0;   // Estimated error for the trigger rate in Hz
} 

// --------------------------------------------------------------------------
//
//  overloaded constructor I
//  fills all member data with initial values and sets the rate of
//  incident showers to ShowRate
//
MHMcRate::MHMcRate (Float_t ShowRate, 
                  const char *name, const char *title)
{

  *fName  = name  ? name  : "MMcTriggerRate";
  *fTitle = title ? title : "Task to calc the collection area ";

  fPartId=0;              // Type of particle
  
  fEnergyMax=0.0;         // Maximum Energy in GeV
  fEnergyMin=1000000.0;   // Minimum Energy in GeV
  
  fThetaMax=0.0;          // Maximum theta angle of run 
  fThetaMin=370.0;        // Minimum theta angle of run
  fPhiMax=0.0;            // Maximum phi angle of run 
  fPhiMin=370.0;          // Minimum phi angle of run
  
  fImpactMax=0.0;         // Maximum impact parameter
  fImpactMin=100000.0;    // Minimum impact parameter
  
  fBackTrig=-1.0;         // Number of triggers from background
  fBackSim=-1.0;          // Number of simulated showers for the background
  
  fSpecIndex=0.0;         // dn/dE = k * e^{- fSpecIndex}
  fFlux0=-1.0;            // dn/dE = fFlux0 * E^{-a}
  
  fShowRate= ShowRate;               // Showers rate in Hz
  fShowRateError=sqrt(ShowRate);     // Estimated error of shower rate in Hz
  
  fTrigRate= -1.0;        // Trigger rate in Hz
  fTrigRateError= -1.0;   // Estimated error for the trigger rate in Hz
} 

// --------------------------------------------------------------------------
//
//  overloaded constructor I
//  fills all member data with initial values and sets the
//  spectral index and the initial flux to SpecIndex and Flux0
//
MHMcRate::MHMcRate (Float_t SpecIndex, Float_t Flux0, 
                  const char *name, const char *title)
{
  *fName  = name  ? name  : "MMcTriggerRate";
  *fTitle = title ? title : "Task to calc the collection area ";

  fPartId=0;              // Type of particle
  
  fEnergyMax=0.0;         // Maximum Energy in GeV
  fEnergyMin=1000000.0;   // Minimum Energy in GeV
  
  fThetaMax=0.0;          // Maximum theta angle of run 
  fThetaMin=370.0;        // Minimum theta angle of run
  fPhiMax=0.0;            // Maximum phi angle of run 
  fPhiMin=370.0;          // Minimum phi angle of run
  
  fImpactMax=0.0;         // Maximum impact parameter
  fImpactMin=100000.0;    // Minimum impact parameter
  
  fBackTrig=-1.0;         // Number of triggers from background
  fBackSim=-1.0;          // Number of simulated showers for the background
  
  fSpecIndex=SpecIndex;   // dn/dE = k * e^{- fSpecIndex}
  fFlux0=Flux0;           // dn/dE = fFlux0 * E^{-a}
  
  fShowRate= -1.0;       
  fShowRateError=0.0;     
  
  fTrigRate= -1.0;        // Trigger rate in Hz
  fTrigRateError= -1.0;   // Estimated error for the trigger rate in Hz
} 

// --------------------------------------------------------------------------
//
//  default constructor
//
MHMcRate::~MHMcRate()
{ 
} 

// --------------------------------------------------------------------------
//
//  set the particle that produces the showers in the athmosphere
//
void MHMcRate:: SetParticle(UShort_t part)
{ 
  fPartId=part;
} 

// --------------------------------------------------------------------------
//
//  Set the information about trigger due only to the background conditions
//
void MHMcRate::SetBackground (Float_t Showers, Float_t Triggers)
{ 
  fBackTrig=Showers;      // Number of triggers from background
  fBackSim=Triggers;      // Number of simulated showers for the background
} 

// --------------------------------------------------------------------------
//
//  update the limits for energy, theta, phi and impact parameter
//
void MHMcRate::Boundaries(Float_t Energy, Float_t Theta, 
                         Float_t Phi, Float_t Impact)
{ 
  // It updates the limit values 

  if(fThetaMax<Theta) fThetaMax=Theta;
  if(fThetaMin>Theta) fThetaMin=Theta;

  if(fPhiMax<Phi) fPhiMax=Phi;
  if(fPhiMin>Phi) fPhiMin=Phi;

  if(fImpactMax<Impact) fImpactMax=Impact;
  if(fImpactMin>Impact) fImpactMin=Impact;

  if(fEnergyMax<Energy) fEnergyMax=Energy;
  if(fEnergyMin>Energy) fEnergyMin=Energy;

} 

// --------------------------------------------------------------------------
//
//  compute the trigger rate
//
void MHMcRate::CalcRate(Float_t trig, Float_t anal, Float_t simu) 
{ 
  // It computes the trigger rate

  // First one computes the rate of incident showers.
  fShowRate=fFlux0/(fSpecIndex-1)*
    (-pow(fEnergyMax,1-fSpecIndex)+pow(fEnergyMin,1-fSpecIndex));

  if (fPartId!=1) 
    fShowRate=fShowRate*(fPhiMax-fPhiMin)*
      (cos(fThetaMax)-cos(fThetaMin));

  fShowRate=fShowRate*PI_NUMBER*(fImpactMax/100.0*fImpactMax/100.0-
                                   fImpactMin/100.0*fImpactMin/100.0);

  fShowRateError=sqrt(fShowRate);

  // Then the trigger rate and its error is evaluated
  if(fBackTrig<0){
    fTrigRateError=sqrt((trig*fShowRate*fShowRate/(simu*simu))+
                    ((1.0-anal*160.0e-9)*
                     (1.0-anal*160.0e-9)*1/
                     (fBackSim*fBackSim*160.0e-9*fBackSim*160.0e-9))); 
    fBackTrig=0;
  }
  else
    fTrigRateError=sqrt((trig*fShowRate*fShowRate/(simu*simu))+
                       ((1.0-anal*160.0e-9)*
                        (1.0-anal*160.0e-9)*fBackTrig/
                        (fBackSim*160.0e-9*fBackSim*160.0e-9))); 
  
  fTrigRate=trig*fShowRate/simu+
      (1.0-anal*160.0e-9)*fBackTrig/
      (fBackSim*160.0e-9);

} 

// --------------------------------------------------------------------------
//
//  print the trigger rate
//
void MHMcRate::Print()
{
  *fLog<<"Incident rate "<<fShowRate<<" Hz "<<endl;
  *fLog<<endl<<"Trigger Rate "<<fTrigRate<<" +- "<<fTrigRateError<<" Hz"<<endl;

}

// --------------------------------------------------------------------------
//
//  draw the trigger rate
//
void MHMcRate::Draw(Option_t* option) 
{ 
  *fLog<<"To be iplemented"<<endl;
}