source: trunk/MagicSoft/Mars/mmontecarlo/MMcEnerThreCalc.cc@ 860

/* ======================================================================== *\
!
! *
! * 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): Javier Lopez 05/2001 (jlopez@ifae.es)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MMcEnerThreCalc
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "MMcEnerThreCalc.h"

#include "MParList.h"

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

#include "MMcEvt.hxx"
#include "MMcTrig.hxx"
#include "MMcEnerThre.h"
#include "MMcEnerHisto.h"

#include <math.h>

ClassImp(MMcEnerThreCalc)

MMcEnerThreCalc::MMcEnerThreCalc (const int dim,
                                  const char* name, const char* title)
{
  *fName  = name  ? name  : "MMcEnerThreCalc";
  *fTitle = title ? title : "Task to calculate the energy threshold from Monte Carlo";

  fDimension = dim;

  if (fDimension > 0)
    fMMcTrig = new MMcTrig * [fDimension];
  else
    fMMcTrig = new MMcTrig * [1];
}

Bool_t MMcEnerThreCalc::PreProcess(MParList* pList)
{
  // connect Monte Carlo data with this task

  char auxname[15]; // string to write container names

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

  if (fDimension == 0)
    {
      fMMcTrig[0] = (MMcTrig*)pList->FindObject("MMcTrig");
      if (!fMMcTrig[0])
        {
          *fLog << dbginf << "MMcTrig not found... aborting." << endl;
          return kFALSE;
        }
    }
  else
    {
      for (int i=0; i<fDimension; i++)
        {
          sprintf(auxname,"MMcTrig;%i.",i+1);
          fMMcTrig[i] = (MMcTrig*)pList->FindObject(auxname);
          if (!fMMcTrig[i])
            {
              *fLog << dbginf << "MMcTrig;" << i+1 << "  not found... aborting." << endl;
              return kFALSE;
            } 
        }
    }

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

Bool_t MMcEnerThreCalc::Process()
{
  const Float_t energy = fMMcEvt->GetEnergy();
  if (fDimension <= 0)
    {
      Int_t trigger[1]; 
      trigger[0] = fMMcTrig[0]->GetFirstLevel();
      if (trigger[0]>0) (*fMMcEnerThre)[0]->Fill(log10(energy),1/energy);
    }
  else
    {
      Int_t *trigger;
      trigger=new Int_t [fDimension];

      for (int i=0; i<fDimension; i++)
        {
          trigger[i] = fMMcTrig[i]->GetFirstLevel() ;
          if (trigger[i]>0) (*fMMcEnerThre)[i]->Fill(log10(energy),1/energy);
        }
    }
  return kTRUE;
}

Bool_t MMcEnerThreCalc::PostProcess()
{
  // fit the energy distribution to get the threshold

  char auxname[15];
  
  if (fDimension <= 0) 
    {
      (*fMMcEnerThre)[0]->Fit("fLogEner0","RQ","",1.,3.);
      (*fMMcEnerThre)[0]->Fit("fLogEner0","RQ","",(*fMMcEnerThre)[0]->GetPeakAtLogE()-2*(*fMMcEnerThre)[0]->GetSigmaAtLogE(),(*fMMcEnerThre)[0]->GetPeakAtLogE()+2*(*fMMcEnerThre)[0]->GetSigmaAtLogE());
      (*fMMcEnerThre)[0]->Fit("fLogEner0","RQ","",(*fMMcEnerThre)[0]->GetPeakAtLogE()-sqrt(2)*(*fMMcEnerThre)[0]->GetSigmaAtLogE(),(*fMMcEnerThre)[0]->GetPeakAtLogE()+sqrt(2)*(*fMMcEnerThre)[0]->GetSigmaAtLogE());
 
    }
  else
    {
      for (int i=0; i<fDimension; i++)
        {
          sprintf (auxname,"fLogEner%i",i);
          
          (*fMMcEnerThre)[i]->Fit(auxname,"RQ","",1.,3.);
          (*fMMcEnerThre)[i]->Fit(auxname,"RQ","",(*fMMcEnerThre)[i]->GetPeakAtLogE()-2*(*fMMcEnerThre)[i]->GetSigmaAtLogE(),(*fMMcEnerThre)[i]->GetPeakAtLogE()+2*(*fMMcEnerThre)[i]->GetSigmaAtLogE());
          (*fMMcEnerThre)[i]->Fit(auxname,"RQ","",(*fMMcEnerThre)[i]->GetPeakAtLogE()-sqrt(2)*(*fMMcEnerThre)[i]->GetSigmaAtLogE(),(*fMMcEnerThre)[i]->GetPeakAtLogE()+sqrt(2)*(*fMMcEnerThre)[i]->GetSigmaAtLogE());
        }
    }
  return kTRUE;
}