source: trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc@ 9301

/* ======================================================================== *\
!
! *
! * This file is part of CheObs, the Modular 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 appears 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,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MSimRandomPhotons
//
//  Simulate poissonian photons. Since the distribution of the arrival time
// differences of these photons is an exonential we can simulate them
// using exponentially distributed time differences between two consecutive
// photons.
//
// FIXME: We should add the wavelength distribution.
//
//  Input Containers:
//   fNameGeomCam [MGeomCam]
//   MPhotonEvent
//   MPhotonStatistics
//   MCorsikaEvtHeader
//   [MCorsikaRunHeader]
//
//  Output Containers:
//   MPhotonEvent
//
//////////////////////////////////////////////////////////////////////////////
#include "MSimRandomPhotons.h"

#include <TRandom.h>

#include "MMath.h"        // RndmExp

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MGeomPix.h"

#include "MPhotonEvent.h"
#include "MPhotonData.h"

#include "MCorsikaRunHeader.h"

ClassImp(MSimRandomPhotons);

using namespace std;

// --------------------------------------------------------------------------
//
//  Default Constructor.
//
MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title)
    : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0),
    fSimulateWavelength(kFALSE),  fNameGeomCam("MGeomCam")
{
    fName  = name  ? name  : "MSimRandomPhotons";
    fTitle = title ? title : "Simulate possonian photons (like NSB or dark current)";
}

// --------------------------------------------------------------------------
//
//  Check for the necessary containers
//
Int_t MSimRandomPhotons::PreProcess(MParList *pList)
{
    fGeom = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam");
    if (!fGeom)
    {
        *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;

        fGeom = (MGeomCam*)pList->FindCreateObj(fNameGeomCam);
        if (!fGeom)
            return kFALSE;
    }

    fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
    if (!fEvt)
    {
        *fLog << err << "MPhotonEvent not found... aborting." << endl;
        return kFALSE;
    }

    fStat = (MPhotonStatistics*)pList->FindObject("MPhotonStatistics");
    if (!fStat)
    {
        *fLog << err << "MPhotonStatistics not found... aborting." << endl;
        return kFALSE;
    }

    /*
    fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader");
    if (!fEvtHeader)
    {
        *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl;
        return kFALSE;
    }*/

    fRunHeader = 0;
    if (fSimulateWavelength)
    {
        fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
        if (!fRunHeader)
        {
            *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
            return kFALSE;
        }
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Check for the necessary containers
//
Int_t MSimRandomPhotons::Process()
{
    // Get array from event container
    // const Int_t num = fEvt->GetNumPhotons();
    //
    // Do not produce pure pedestal events!
    // if (num==0)
    //    return kTRUE;

    // Get array from event container
    // FIXME: Use statistics container instead
    const UInt_t npix = fGeom->GetNumPixels();

    // This is the possible window in which the triggered digitization
    // may take place.
    const Double_t start = fStat->GetTimeFirst();
    const Double_t end   = fStat->GetTimeLast();

    // Loop over all pixels
    for (UInt_t idx=0; idx<npix; idx++)
    {
        // Scale the rate with the pixel size. The rate must
        // always be given for the pixel with index 0.
        const Double_t avglen = 1./(fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()/100.);

        // Start producing photons at time "start"
        Double_t t = start;
        while (1)
        {
            // Get a random time for the photon.
            // The differences are exponentially distributed.
            t += MMath::RndmExp(avglen);

            // Check if we reached the end of the useful time window
            if (t>end)
                break;

            // Add a new photon
            // FIXME: SLOW!
            MPhotonData &ph = fEvt->Add(); 

            // Set source to NightSky, time to t and tag to pixel index
            ph.SetPrimary(MMcEvtBasic::kNightSky);
            ph.SetWeight();
            ph.SetTime(t);
            ph.SetTag(idx);

            // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant)  FIXME: Reset?

            if (fRunHeader)
            {
                const Float_t wmin = fRunHeader->GetWavelengthMin();
                const Float_t wmax = fRunHeader->GetWavelengthMax();

                ph.SetWavelength(TMath::Nint(gRandom->Uniform(wmin, wmax)));
            }
        }
    }

    // Re-sort the photons by time!
    fEvt->Sort(kTRUE);

    // Update maximum index
    fStat->SetMaxIndex(npix-1);

    // Shrink
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read the parameters from the resource file.
//
//    FrequencyFixed: 40
//    FrequencyNSB:   40
//
// The frequency is given in units fitting the units of the time.
// Usually the time is given in nanoseconds thus, e.g., 40 means 40MHz
//
Int_t MSimRandomPhotons::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "FrequencyFixed", print))
    {
        rc = kTRUE;
        fFreqFixed = GetEnvValue(env, prefix, "FrequencyFixed", fFreqFixed);
    }

    if (IsEnvDefined(env, prefix, "FrequencyNSB", print))
    {
        rc = kTRUE;
        fFreqNSB = GetEnvValue(env, prefix, "FrequencyNSB", fFreqNSB);
    }

    return rc;
}