Index: trunk/MagicSoft/Mars/msimcamera/MSimBundlePhotons.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimBundlePhotons.cc (revision 9241) +++ trunk/MagicSoft/Mars/msimcamera/MSimBundlePhotons.cc (revision 9241) @@ -0,0 +1,203 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz 1/2009 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSimBundlePhotons +// +// This task sets a new tag (index) for all photons in a list MPhotonEvent +// based on a look-up table. +// +// Input: +// MPhotonEvent +// MPhotonStatistics +// +// Output +// MPhotonEvent +// +////////////////////////////////////////////////////////////////////////////// +#include "MSimBundlePhotons.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MArrayI.h" + +#include "MParList.h" + +#include "MPhotonEvent.h" +#include "MPhotonData.h" + +ClassImp(MSimBundlePhotons); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MSimBundlePhotons::MSimBundlePhotons(const char* name, const char *title) +: fEvt(0), fStat(0)//, fFileName("mreflector/dwarf-apdmap.txt") +{ + fName = name ? name : "MSimBundlePhotons"; + fTitle = title ? title : "Task to bundle (re-index) photons according to a look-up table"; +} + +// -------------------------------------------------------------------------- +// +// Check for the needed containers and read the oruting table. +// +Int_t MSimBundlePhotons::PreProcess(MParList *pList) +{ + 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; + } + + if (fFileName.IsNull()) + return kSKIP; + + // Read the look-up table + if (fLut.ReadFile(fFileName)<0) + return kFALSE; + + // If the table is empty remove this task from the tasklist + if (fLut.IsEmpty()) + return kSKIP; + + // Now invert the tablee. Otherwise we have to do a lot of + // searching for every index. + fLut.Invert(); + + // Make sure that each line has exactly one row + if (!fLut.HasConstantLength() && fLut.GetMaxEntries()!=1) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Re-index all photons according to the look-up table. +// +Int_t MSimBundlePhotons::Process() +{ + // Make sure that we don't get a seg-fault + /* + if (fStat->GetMaxIndex()>=fLut.GetEntriesFast()) + { + *fLog << err; + *fLog << "ERROR - MSimBundlePhotons::Process: Maximum pixel index stored" << endl; + *fLog << " in tag of MPhotonData exceeds the look-up table length." << endl; + return kERROR; + }*/ + + // FIXME: Add a range check comparing the min and max tag with + // the number of entries in the routing table + + // Get total number of photons + const Int_t num = fEvt->GetNumPhotons(); + + // If there are no photons we can do nothing + if (num==0) + return kTRUE; + + // Get maximum index allowed + const Int_t max = fLut.GetEntriesFast(); + + // Initialize a counter for the final number of photons. + Int_t cnt=0; + + // Loop over all photons + for (Int_t i=0; i=max) + continue; + + // Get the routing assigned to this index + const MArrayI &row = fLut.GetRow(tag); + + // Sanity check: Check if we were routed to a + // valid entry if not throw away this photon. + if (row.GetSize()==0) + continue; + + // Get corresponding entry from routing table + const Int_t &idx = row[0]; + + // Check if we were routed to a valid entry + // if not throw away this photon. + if (idx<0) + continue; + + // Set Tag to new index + ph.SetTag(idx); + + // Copy photon to its now position in array and increade counter + (*fEvt)[cnt++] = ph; + } + + // Shrink the list of photons to its new size + fEvt->Shrink(cnt); + + // Set new maximum index (Note, that this is the maximum index + // available in the LUT, which does not necessarily correspond + // to, e.g., the number of pixels although it should) + fStat->SetMaxIndex(fLut.GetMaxIndex()); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// FileName: lut.txt +// +Int_t MSimBundlePhotons::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "FileName", print)) + { + rc = kTRUE; + fFileName = GetEnvValue(env, prefix, "FileName", fFileName); + } + + return rc; +} Index: trunk/MagicSoft/Mars/msimcamera/MSimBundlePhotons.h =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimBundlePhotons.h (revision 9241) +++ trunk/MagicSoft/Mars/msimcamera/MSimBundlePhotons.h (revision 9241) @@ -0,0 +1,40 @@ +#ifndef MARS_MSimBundlePhotons +#define MARS_MSimBundlePhotons + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MLut +#include "MLut.h" +#endif + +class MParList; +class MPhotonEvent; +class MPhotonStatistics; + +class MSimBundlePhotons: public MTask +{ +private: + MPhotonEvent *fEvt; //! Event storing the photons + MPhotonStatistics *fStat; //! Event statistics needed for crosschecks + + TString fFileName; // File to from which to read the lut + MLut fLut; // Look-up table + + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + + // MTask + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MSimBundlePhotons(const char *name=NULL, const char *title=NULL); + + void SetFileName(const char *name) { fFileName = name; } + + ClassDef(MSimBundlePhotons, 0) // Task to bundle (re-index) photons according to a look-up table +}; + +#endif Index: trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc (revision 9241) +++ trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc (revision 9241) @@ -0,0 +1,232 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz 12/2008 +! +! Copyright: 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 + +#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; idxend) + 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.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; +} Index: trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h (revision 9241) +++ trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h (revision 9241) @@ -0,0 +1,50 @@ +#ifndef MARS_MSimRandomPhotons +#define MARS_MSimRandomPhotons + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MGeomCam; +class MParList; +class MPhotonEvent; +class MPhotonStatistics; +//class MCorsikaEvtHeader; +class MCorsikaRunHeader; + +class MSimRandomPhotons : public MTask +{ +private: + MGeomCam *fGeom; //! container with the geometry + MPhotonEvent *fEvt; //! Event storing the photons + MPhotonStatistics *fStat; //! Container storing evenet statistics +// MCorsikaEvtHeader *fEvtHeader; //! Header storing event information + MCorsikaRunHeader *fRunHeader; //! Header storing run information + + + // FIXME: Make this a single number per Pixel/APD + Double_t fFreqFixed; // [1/ns] A fixed frequency per pixel + Double_t fFreqNSB; // [1/ns/cm^2] A frequency depending on area + + Bool_t fSimulateWavelength; + + TString fNameGeomCam; + + // MTask + Int_t PreProcess(MParList *pList); + Int_t Process(); + + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + +public: + MSimRandomPhotons(const char *name=NULL, const char *title=NULL); + + void SetFreq(Float_t fnsb, Float_t fdc) { fFreqNSB=fnsb; fFreqFixed=fdc; } + + void SetNameGeomCam(const char *name="MGeomCam") { fNameGeomCam = name; } + + ClassDef(MSimRandomPhotons, 0) // Simulate possonian photons (like NSB or dark current) +}; + +#endif