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

/* ======================================================================== *\
!
! *
! * 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
//   AccidentalPhotonRate [MPedestalCam]
//
//////////////////////////////////////////////////////////////////////////////
#include "MSimRandomPhotons.h"

#include <TRandom.h>

#include "MMath.h"        // RndmExp

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MGeom.h"

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

#include "MPedestalCam.h"
#include "MPedestalPix.h"

#include "MCorsikaRunHeader.h"

#include "MSpline3.h"
#include "MParSpline.h"
#include "MReflector.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),
    fRates(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->FindObject("MGeomCam");
        if (!fGeom)
        {
            *fLog << err << "MGeomCam not found... aborting." << endl;
            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;
    }

    fRates = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "AccidentalPhotonRates");
    if (!fRates)
        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;
        }
    }

    MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
    if (!r)
    {
        *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
        return kFALSE;
    }

    const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline");
    const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesAngularAcceptance",    "MParSpline");
    const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity",        "MParSpline");

    const Double_t d2  = fGeom->GetCameraDist()*fGeom->GetCameraDist();
    const Double_t pde = s1 && s1->GetSpline() ? s1->GetSpline()->Integral() : 1;
    const Double_t sr  = s2 && s2->GetSpline() ? s2->GetSpline()->IntegralSolidAngle() : 1;
    const Double_t mir = s3 && s3->GetSpline() ? s3->GetSpline()->Integral() : 1;
    const Double_t Ar  = r->GetA()/1e4;

    // Conversion factor to convert pixel area to steradians (because it
    // is a rather small area we can assume it is flat)
    const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad();

    // Multiply all relevant efficiencies
    MParSpline *s4 = (MParSpline*)s1->Clone();
    s4->Multiply(*s3->GetSpline());

    const Double_t nm = s4 && s4->GetSpline() ? s4->GetSpline()->Integral() : 1;

    delete s4;

    // /100 to convert the pixel area from mm^2 to cm^2
    fScale =  nm * TMath::Min(Ar, sr*d2) * conv*conv;

    *fLog << inf;
    *fLog << "Effective cone acceptance:      " << Form("%.2f", sr*d2) << "m^2" << endl;
    *fLog << "Reflector area:                 " << Form("%.2f", Ar) << "m^2" << endl;
    *fLog << "Resulting eff. collection area: " << Form("%.1f", TMath::Min(Ar, sr*d2)) << "m^2" << endl;
    *fLog << "Eff. wavelength band (PDE):     " << Form("%.1f", pde) << "nm" << endl;
    *fLog << "Eff. wavelength band (Mirror):  " << Form("%.1f", mir) << "nm" << endl;
    *fLog << "Eff. wavelength band (PDE+MIR): " << Form("%.1f", nm) << "nm" << endl;
    *fLog << "Pixel area of " << fNameGeomCam << "[0]: " << Form("%.2e", (*fGeom)[0].GetA()*conv*conv) << "sr" << endl;
    //*fLog << "Effective angular acceptance:     " << sr << "sr" << endl;
    //*fLog << "Resulting NSB frequency:           " << fFreqNSB*nm*Ar*1000 << "MHz/sr" << endl;
    *fLog << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", fFreqNSB*(*fGeom)[0].GetA()*fScale*1000) << "MHz" << endl;

    // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
    // Set NumPheFromDNSB

    // # Number of photons from the diffuse NSB (nphe / ns 0.1*0.1 deg^2 239 m^2) and
    // nsb_mean 0.20
    // Magic pixel: 0.00885361 deg
    // dnsbpix = 0.2*50/15
    // ampl = MMcFadcHeader->GetAmplitud()
    // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio)

    return kTRUE;
}

Bool_t MSimRandomPhotons::ReInit(MParList *pList)
{
    // Overwrite the default set by MGeomApply
    fRates->Init(*fGeom);
    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.
        const Double_t rate = fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()*fScale;

        (*fRates)[idx].SetPedestal(rate);

        // Calculate the average distance between two consequtive photons
        const Double_t avglen = 1./rate;

        // 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: 0.040
//    FrequencyNSB:   0.040
//
// The fixed frequency is given in units fitting the units of the time.
// Usually the time is given in nanoseconds thus, e.g., 0.040 means 40MHz.
//
// The FrequencyNSB is scaled by the area of the pixel in cm^2. Therefore
// 0.040 would mean 40MHz/cm^2
//
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;
}