Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 9241) +++ trunk/MagicSoft/Mars/Changelog (revision 9242) @@ -18,4 +18,20 @@ -*-*- END OF LINE -*-*- + + + 2009/01/23 Thomas Bretz + + * msimcamera/MSimAPD.[h,cc], msimcamera/MSimBundlePhotons.[h,cc], + msimcamera/MSimRandomPhotons.[h,cc]: + - added + + * melectronics/Makefile, melectronics/MAnalogChannels.[h,cc], + melectronics/ElectronicsIncl.h, melectronics/ElectronicsLinkDef.h, + melectronics/MDigitalSignal.[h,cc], + melectronics/MAnalogSignal.[h,cc], + melectronics/MAvalanchePhotoDiode.[h,cc]: + - added + + 2009/01/22 Thomas Bretz Index: trunk/MagicSoft/Mars/melectronics/ElectronicsIncl.h =================================================================== --- trunk/MagicSoft/Mars/melectronics/ElectronicsIncl.h (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/ElectronicsIncl.h (revision 9242) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: trunk/MagicSoft/Mars/melectronics/ElectronicsLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/melectronics/ElectronicsLinkDef.h (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/ElectronicsLinkDef.h (revision 9242) @@ -0,0 +1,12 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class APD+; +#pragma link C++ class MAnalogSignal+; +#pragma link C++ class MDigitalSignal+; +#pragma link C++ class MAnalogChannels+; + +#endif Index: trunk/MagicSoft/Mars/melectronics/MAnalogChannels.cc =================================================================== --- trunk/MagicSoft/Mars/melectronics/MAnalogChannels.cc (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MAnalogChannels.cc (revision 9242) @@ -0,0 +1,148 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MAnalogChannels +// +// A collection of MAnalogSignals with a simplified access interface. +// +////////////////////////////////////////////////////////////////////////////// +#include "MAnalogChannels.h" + +#include + +#include "MAnalogSignal.h" + +ClassImp(MAnalogChannels); + +using namespace std; + +// ------------------------------------------------------------------------ +// +// Default constructor +// +MAnalogChannels::MAnalogChannels(const char *name, const char *title) +: fArray(0), fValidRangeMin(0), fValidRangeMax(0) +{ + fName = name ? name : "MAnalogChannels"; + fTitle = title ? title : "Parameter container for a collection of analog signals"; +} + +// ------------------------------------------------------------------------ +// +// Constructor. Calls Init(n, len) +// +MAnalogChannels::MAnalogChannels(Int_t n, Int_t len, const char *name, const char *title) +: fArray(0), fValidRangeMin(0), fValidRangeMax(0) +{ + fName = name ? name : "MAnalogChannels"; + fTitle = title ? title : "Parameter container for a collection of analog signals"; + + Init(n, len); +} + +// ------------------------------------------------------------------------ +// +// delete fArray +// +void MAnalogChannels::Clear(Option_t *o) +{ + if (fArray) + { + delete fArray; + fArray = 0; + } +} + +// ------------------------------------------------------------------------ +// +// Initializes n channels with a length of n sampling points. +// +// FIXME: We could improve speed by only increasing the arrays? +// +void MAnalogChannels::Init(UInt_t n, UInt_t len) +{ + Clear(); + + fArray = new TObjArray(n); + fArray->SetOwner(); + + for (UInt_t i=0; iAddAt(new MAnalogSignal(len), i); +} + +// ------------------------------------------------------------------------ +// +// return an unchecked reference to the i-th channel. +// +MAnalogSignal &MAnalogChannels::operator[](UInt_t i) +{ + return *static_cast(fArray->UncheckedAt(i)); +} + +// ------------------------------------------------------------------------ +// +// return an unchecked pointer to the i-th channel. +// +MAnalogSignal *MAnalogChannels::operator()(UInt_t i) +{ + return static_cast(fArray->UncheckedAt(i)); +} + +// ------------------------------------------------------------------------ +// +// Return an unchecked const reference to the i-th channel. +// +const MAnalogSignal &MAnalogChannels::operator[](UInt_t i) const +{ + return *static_cast(fArray->UncheckedAt(i)); +} + +// ------------------------------------------------------------------------ +// +// Return an unchecked const pointer to the i-th channel. +// +const MAnalogSignal *MAnalogChannels::operator()(UInt_t i) const +{ + return static_cast(fArray->UncheckedAt(i)); +} + +// ------------------------------------------------------------------------ +// +// Return the number of channels. +// +Int_t MAnalogChannels::GetNumChannels() const +{ + return fArray ? fArray->GetEntriesFast() : 0; +} + +// ------------------------------------------------------------------------ +// +// Return the number of samples per channel +// +Int_t MAnalogChannels::GetNumSamples() const +{ + return operator()(0) ? operator()(0)->GetSize() : 0; +} Index: trunk/MagicSoft/Mars/melectronics/MAnalogChannels.h =================================================================== --- trunk/MagicSoft/Mars/melectronics/MAnalogChannels.h (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MAnalogChannels.h (revision 9242) @@ -0,0 +1,48 @@ +#ifndef MARS_MAnalogChannels +#define MARS_MAnalogChannels + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class TObjArray; +class MAnalogSignal; + +class MAnalogChannels : public MParContainer +{ +private: + TObjArray *fArray; + + Int_t fValidRangeMin; + Int_t fValidRangeMax; + +public: + MAnalogChannels(const char *name=0, const char *title=0); + MAnalogChannels(Int_t n, Int_t len, const char *name=0, const char *title=0); + ~MAnalogChannels() + { + Clear(); + } + + void Clear(Option_t *o=""); + + void Init(UInt_t n, UInt_t len); + + MAnalogSignal &operator[](UInt_t i); + MAnalogSignal *operator()(UInt_t i); + + const MAnalogSignal &operator[](UInt_t i) const; + const MAnalogSignal *operator()(UInt_t i) const; + + Int_t GetNumChannels() const; + Int_t GetNumSamples() const; + + void SetValidRange(Int_t min, Int_t max) { fValidRangeMin=min; fValidRangeMax=max; } + + Int_t GetValidRangeMin() const { return fValidRangeMin; } + Int_t GetValidRangeMax() const { return fValidRangeMax; } + + ClassDef(MAnalogChannels, 1) // Parameter container for a collection of analog signals +}; + +#endif Index: trunk/MagicSoft/Mars/melectronics/MAnalogSignal.cc =================================================================== --- trunk/MagicSoft/Mars/melectronics/MAnalogSignal.cc (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MAnalogSignal.cc (revision 9242) @@ -0,0 +1,233 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MAnalogSignal +// +// This is the equivalent to an analog signal. The signal is stored by +// a sampling in equidistant bins. +// +////////////////////////////////////////////////////////////////////////////// +#include "MAnalogSignal.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MSpline3.h" +#include "MDigitalSignal.h" + +#include "MExtralgoSpline.h" + +ClassImp(MAnalogSignal); + +using namespace std; + +// ------------------------------------------------------------------------ +// +// Set the array length and the length of the buffers. +// +void MAnalogSignal::Set(UInt_t n) +{ + // FIXME: Maybe we move this before initializing the spline + // with a check? + fDer1.Set(n); + fDer2.Set(n); + + MArrayF::Set(n); +} + +// ------------------------------------------------------------------------ +// +// Evaluate the spline an add the result between t+xmin and t+xmax +// (xmin and xmax are the limits of the spline) to the signal. +// The spline is evaluated at the bin-center of the analog signal +// and multiplied by f. +// +void MAnalogSignal::AddPulse(const MSpline3 &spline, Float_t t, Float_t f) +{ + // FIXME: This could be improved using a MExtralgoSpline with + // the identical stepping as the signal and we could use + // the integral instead of the pure y-value if we want. + + // Both in units of the sampling frequency + const Float_t start = t+spline.GetXmin(); + const Float_t end = t+spline.GetXmax(); + + /*const*/ Int_t first = TMath::CeilNint(start); + /*const*/ UInt_t last = TMath::CeilNint(end); // Ceil:< Floor:<= + + if (first<0 || last>GetSize()) + { + gLog << err << "ERROR - AddPulse: Out of bounds, "; + gLog << "Win=[" << first << "," << last << "] N=" << GetSize() << " t=" << t; + gLog << " Spline=[" << spline.GetXmin() << "," << spline.GetXmax() << "]" << endl; + + if (first<0) + first=0; + if (last>GetSize()) + last=GetSize(); + } + + Float_t *arr = GetArray(); + for (UInt_t i=first; iExp(num); + AddPulse(spline, d); + } +} + +// ------------------------------------------------------------------------ +// +// Add a random gaussian with amplitude and offset to every bin +// of the analog signal. The default offset is 0. The default amplitude 1. +// +void MAnalogSignal::AddGaussianNoise(Float_t amplitude, Float_t offset) +{ + for (Float_t *ptr = GetArray(); ptrGaus(offset, amplitude); +} + +// ------------------------------------------------------------------------ +// +// The signal is evaluated using the spline MExtralgoSpline. +// Searching upwards from the beginning all points are calculated at +// which the spline is equal to threshold. After a rising edge +// a leading edge is searched. From this an MDigitalSignal is +// created and added to an newly created TObjArray. If len<0 then +// the signal length is equal to the time above threshold, otherwise +// the length is fixed to len. The start of the digital signal is the +// rising edge. If due to fixed length two digital signal overlap the +// digital signals are combined into one signal. +// +// For numerical reasons we have to avoid to find the same x-value twice. +// Therefor a "dead-time" of 1e-4 is implemented after each edge. +// +// The user is responsible of deleting the TObjArray. +// +TObjArray *MAnalogSignal::Discriminate(Float_t threshold, Float_t len) const +{ + TObjArray *ttl = new TObjArray; + ttl->SetOwner(); + + // The time after which we start searching for a falling or leading + // edge at threshold after a leading or falling edge respectively. + // This value has mainly numerical reasons. If starting the search + // too early we might end up in an endless loop finding the same + // value again and again. This just means that a glitch above or + // below the threshold which is shorter than this value can + // stay unnoticed. This shouldn't hurt. + const Double_t deadtime = 1e-4; + + // FIXME: Are local maximum/minima handled correctly? + + const MExtralgoSpline sp(GetArray(), fN, fDer1.GetArray(), fDer2.GetArray()); + + Double_t x1 = 0; + Double_t x2 = 0; + + while (1) + { + x1 = sp.SearchYup(x2+deadtime, threshold); + if (x1<0) + break; + + const Bool_t rising = sp.Deriv1(x1)>0; + if (!rising) + { + // The last value might just have been a local max/min + // or its period was short than 1e-4 + x2 = x1; + //gLog << warn << "Rising edge expected at " << x1 << " (after " << x2 << ", N=" << fN << ")" << endl; + continue; + } + + x2 = sp.SearchYup(x1+deadtime, threshold); + if (x2<0) + break; + + const Bool_t falling = sp.Deriv1(x2)<=0; + if (!falling) + { + // The last value might just have been a local max/min + // or its period was short than 1e-4 + x1 = x2; + //gLog << warn << "Falling edge expected at " << x2 << " (after " << x1 << " N=" << fN << ")" << endl; + continue; + } + + MDigitalSignal *sig = new MDigitalSignal(x1, len>0?len:x2-x1); + + // In case of a fixed length we have to check for possible overlapping + if (len>0 && ttl->GetEntriesFast()>0) + { + // FIXME: What if in such a case the electronics is just dead? + MDigitalSignal *last = static_cast(ttl->Last()); + // Combine both signals to one if they overlap + if (last->Combine(*sig)) + { + // Both signals overlap and have been combined into the existing one + delete sig; + continue; + } + // The signals don't overlap we add the new signal as usual + } + + // Add the new signal to the list of signals + ttl->Add(sig); + } + + if (x1>=0) + gLog << warn << "Falling edge expected before end at N=" << fN << "!" << endl; + + return ttl; +} Index: trunk/MagicSoft/Mars/melectronics/MAnalogSignal.h =================================================================== --- trunk/MagicSoft/Mars/melectronics/MAnalogSignal.h (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MAnalogSignal.h (revision 9242) @@ -0,0 +1,33 @@ +#ifndef MARS_MAnalogSignal +#define MARS_MAnalogSignal + +#ifndef MARS_MArrayF +#include "MArrayF.h" +#endif + +class MSpline3; + +class MAnalogSignal : public MArrayF/*TObject*/ +{ +private: + MArrayF fDer1; //! Buffer for the derivatives of the corresponding spline + MArrayF fDer2; //! Buffer for the derivatives of the corresponding spline + +public: + MAnalogSignal(UInt_t n) { Set(n); } + + void Set(UInt_t n); + void AddPulse(const MSpline3 &spline, Float_t t, Float_t f=1); + void AddSignal(const MAnalogSignal &s); + + // Deprecated. Use MSimRandomPhotons instead + void AddRandomPulses(const MSpline3 &spline, Float_t num); + + void AddGaussianNoise(Float_t amplitude=1, Float_t offset=0); + + TObjArray *Discriminate(Float_t threshold, Float_t len=-1) const; + + ClassDef(MAnalogSignal, 1) // Storage class for an analog signal +}; + +#endif Index: trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc =================================================================== --- trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc (revision 9242) @@ -0,0 +1,304 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// APD +// +// All times in this class are relative times. Therefor the unit for the +// time is not intrinsically fixed. In fact the dead-time and recovery- +// time given in the constructor must have the same units. This is what +// defines the unit of the times given in the function and the unit of +// rates given. +// For example, if recovery and dead time are given in nanoseconds the +// all times must be in nanoseconds and rates are given per nanosecond, +// i.e. GHz. +// +// Hamamatsu 30x30 cells: APD(30, 0.2, 3, 35) +// Hamamatsu 60x60 cells: APD(60, 0.2, 3, 8.75) +// +////////////////////////////////////////////////////////////////////////////// +#include "MAvalanchePhotoDiode.h" + +#include + +#include "MMath.h" + +ClassImp(APD); + +using namespace std; + +/* +class MyProfile : public TProfile2D +{ +public: + void AddBinEntry(Int_t cell) { fBinEntries.fArray[cell]++; } +}; +*/ + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +// n is the number od cells in x or y. The APD is assumed to +// be square. +// prob is the crosstalk probability, i.e., the probability that a +// photon which produced an avalanche will create another +// photon in a neighboring cell +// dt is the deadtime, i.e., the time in which the APD cell will show +// no response to a photon after a hit +// rt is the recovering tims, i.e. the exponential (e^(-dt/rt)) +// with which the cell is recovering after being dead +// +// prob, dt and ar can be set to 0 to switch the effect off. +// 0 is also the dfeault for all three. +// +APD::APD(Int_t n, Float_t prob, Float_t dt, Float_t rt) + : fHist("APD", "", n, 0.5, n+0.5, n, 0.5, n+0.5), + fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt) +{ + fHist.SetDirectory(0); +} + +// -------------------------------------------------------------------------- +// +// This is the recursive implementation of a hit. If a photon hits a cell +// at x and y (must be a valid cell!) at time t, at first we check if the +// cell is still dead. If it is not dead we calculate the signal height +// from the recovery time. Now we check with the crosstalk probability +// whether another photon is created. If another photon is created we +// calculate randomly which of the four neighbor cells are hit. +// If the cell is outside the APD the photon is ignored. As many +// new photons are created until our random number is below the crosstak- +// probability. +// +// The total height of the signal (in units of photons) is returned. +// Note, that this can be a fractional number. +// +// This function looks a bit fancy accessing the histogram and works around +// a few histogram functions. This is a speed optimization which works +// around a lot of sanity checks which are obsolete in our case. +// +// The default time is 0. +// +Float_t APD::HitCellImp(Int_t x, Int_t y, Float_t t) +{ + // if (x<1 || x>fHist.GetNbinsX() || + // y<1 || y>fHist.GetNbinsY()) + // return 0; + + // const Int_t cell = fHist.GetBin(x, y); + const Int_t cell = x + (fHist.GetNbinsX()+2)*y; + + // This is the fastes way to access the bin-contents in fArray + Float_t &cont = fHist.GetArray()[cell]; + + // const Double_t dt = t-fHist.GetBinContent(x, y)-fDeadTime; // + const Float_t dt = t-cont-fDeadTime; + + // Photons within the dead time are just ignored + if (/*hx.GetBinContent(x,y)>0 &&*/ dt<=0) + return 0; + + // Signal height (in units of one photon) produced after dead time + const Float_t weight = fRecoveryTime<=0 ? 1 : 1.-exp(-dt/fRecoveryTime); + + cont = t; // fHist.SetBinContent(x, y, t) + + // Counter for the numbers of produced photons + Float_t n = weight; + + /* + // Check if a photon in a neighboring cell is produced (crosstalk) + while (gRandom->Rndm()Integer(4)) + { + case 0: x++; if (x>fHist.GetNbinsX()) continue; break; + case 1: x--; if (x<1) continue; break; + case 2: y++; if (y>fHist.GetNbinsY()) continue; break; + case 3: y--; if (y<1) continue; break; + } + + n += HitCellImp(x, y, t); + } + */ + + //for (int i=0; i<1; i++) + while (1) + { + const Double_t rndm = gRandom->Rndm(); + if (rndm>=fCrosstalkProb) + break; + + // We can re-use the random number becuase it is uniformely + // distributed. This saves cpu power + const Int_t dir = TMath::FloorNint(4*rndm/fCrosstalkProb); + + // Get a random neighbor which is hit. + switch (dir) + { + case 0: if (x1) n += HitCellImp(x-1, y, t); break; + case 2: if (y1) n += HitCellImp(x, y-1, t); break; + /* + case 0: x++; if (x>fHist.GetNbinsX()) continue; break; + case 1: x--; if (x<1) continue; break; + case 2: y++; if (y>fHist.GetNbinsY()) continue; break; + case 3: y--; if (y<1) continue; break;*/ + } + + // In the unlikely case the calculated direction is out-of-range, + // i.e. <0 or >3, we would just try to fill the same cell again which + + //n += HitCellImp(x, y, t); + } + + return n; +} + +// -------------------------------------------------------------------------- +// +// Check if x and y is a valid cell. If not return 0, otherwise +// HitCelImp(x, y, t) +// +// The default time is 0. +// +Float_t APD::HitCell(Int_t x, Int_t y, Float_t t) +{ + if (x<1 || x>fHist.GetNbinsX() || + y<1 || y>fHist.GetNbinsY()) + return 0; + + return HitCellImp(x, y, t); +} + +// -------------------------------------------------------------------------- +// +// Determine randomly (uniformly) a cell which was hit. Return +// HitCellImp for this cell and the given time. +// +// The default time is 0. +// +Float_t APD::HitRandomCell(Float_t t) +{ + const UInt_t nx = fHist.GetNbinsX(); + const UInt_t ny = fHist.GetNbinsY(); + + const UInt_t idx = gRandom->Integer(nx*ny); + + const UInt_t x = idx%nx; + const UInt_t y = idx/nx; + + return HitCellImp(x+1, y+1, t); +} + +// -------------------------------------------------------------------------- +// +// Sets all cells with a contents whihc is well before the time t such that +// the chip is "virgin". Therefore all cells are set to a time which +// is twice the deadtime before the given time and 1000 times the recovery +// time. +// +// If deadtime and recovery time are 0 then t-1 is set. +// +// The default time is 0. +// +void APD::FillEmpty(Float_t t) +{ + const Int_t n = (fHist.GetNbinsX()+2)*(fHist.GetNbinsY()+2); + + const Double_t tm = fDeadTime<=0 && fRecoveryTime<=0 ? t-1 : t-2*fDeadTime-1000*fRecoveryTime; + + for (int i=0; iExp(n/rate) with n being +// the total number of cells. +// The default time is 0. +// +void APD::FillRandom(Float_t rate, Float_t t) +{ + FillEmpty(t); + + const Int_t nx = fHist.GetNbinsX(); + const Int_t ny = fHist.GetNbinsY(); + + const Double_t f = (nx*ny)/rate; + + // FIXME: This is not perfect, is it? What about the dead time? + + for (int x=1; x<=nx; x++) + for (int y=1; y<=ny; y++) + { + HitCellImp(x, y, t-MMath::RndmExp(f)); + } +} + +// -------------------------------------------------------------------------- +// +// Retunrs the number of cells which have a time t<=fDeadTime, i.e. which are +// dead. +// The default time is 0. +// +Int_t APD::CountDeadCells(Float_t t) const +{ + const Int_t nx = fHist.GetNbinsX(); + const Int_t ny = fHist.GetNbinsY(); + + Int_t n=0; + for (int x=1; x<=nx; x++) + for (int y=1; y<=ny; y++) + if ((t-fHist.GetBinContent(x, y))<=fDeadTime) + n++; + + return n; +} + +// -------------------------------------------------------------------------- +// +// Returs the number of cells which have a time t<=fDeadTime+fRecoveryTime. +// The default time is 0. +// +Int_t APD::CountRecoveringCells(Float_t t) const +{ + const Int_t nx = fHist.GetNbinsX(); + const Int_t ny = fHist.GetNbinsY(); + + Int_t n=0; + for (int x=1; x<=nx; x++) + for (int y=1; y<=ny; y++) + { + Float_t dt = t-fHist.GetBinContent(x, y); + if (dt>fDeadTime && dt<=fDeadTime+fRecoveryTime) + n++; + } + return n; +} Index: trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.h =================================================================== --- trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.h (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.h (revision 9242) @@ -0,0 +1,34 @@ +#ifndef MARS_MAvalanchePhotoDiode +#define MARS_MAvalanchePhotoDiode + +#ifndef ROOT_TH2 +#include +#endif + +class APD : public TObject +{ +private: + TH2F fHist; + + Float_t fCrosstalkProb; // Probability that a converted photon creates another one in a neighboring cell + Float_t fDeadTime; // Deadtime of a single cell after a hit + Float_t fRecoveryTime; // Recoverytime after Deadtime (1-exp(-t/fRecoveryTime) + + Float_t HitCellImp(Int_t x, Int_t y, Float_t t=0); + +public: + APD(Int_t n, Float_t prob=0, Float_t dt=0, Float_t rt=0); + + Float_t HitCell(Int_t x, Int_t y, Float_t t=0); + Float_t HitRandomCell(Float_t t=0); + + void FillEmpty(Float_t t=0); + void FillRandom(Float_t rate, Float_t t=0); + + Int_t CountDeadCells(Float_t t=0) const; + Int_t CountRecoveringCells(Float_t t=0) const; + + ClassDef(APD, 1) // An object representing a Geigermode APD +}; + +#endif Index: trunk/MagicSoft/Mars/melectronics/MDigitalSignal.cc =================================================================== --- trunk/MagicSoft/Mars/melectronics/MDigitalSignal.cc (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MDigitalSignal.cc (revision 9242) @@ -0,0 +1,129 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MDigitalSignal +// +// A digital signal with a start time and length. +// +////////////////////////////////////////////////////////////////////////////// +#include "MDigitalSignal.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MDigitalSignal); + +using namespace std; + +// ------------------------------------------------------------------------ +// +// Initialize a new signal with a coincidence (and) of the two given +// signals. +// +MDigitalSignal::MDigitalSignal(const MDigitalSignal &ttl1, const MDigitalSignal &ttl2) +{ + const Double_t new0 = TMath::Max(ttl1.fStart, ttl2.fStart); + const Double_t new1 = TMath::Min(ttl1.fStart+ttl1.fLength, ttl2.fStart+ttl2.fLength); + + fStart = new0; + fLength = new1-new0; +} + +// ------------------------------------------------------------------------ +// +// Compare the start time of two signal. +// +// Returns: +// 0) Both signals have identical start time +// 1) the argument starts earlier than this +// -1) The argument starts later than this +// +Int_t MDigitalSignal::Compare(const TObject *obj) const +{ + const MDigitalSignal &ttl = *static_cast(obj); + + if (ttl.fStart>fStart) + return -1; + + if (ttl.fStart(obj); + + return fStart<=ttl.fStart+ttl.fLength && ttl.fStart<=fStart+fLength; +} + +// ------------------------------------------------------------------------ +// +// Combine a new (overlapping) signal with this signal to one signal. +// If they don't overlap kFALSE is returned, kTRUE otherwise. +// +Bool_t MDigitalSignal::Combine(const TObject &obj) +{ + if (!Overlap(obj)) + return kFALSE; + + const MDigitalSignal &ttl = static_cast(obj); + + const Double_t new0 = TMath::Min(fStart, ttl.fStart); + const Double_t new1 = TMath::Max(fStart+fLength, ttl.fStart+ttl.fLength); + + fStart = new0; + fLength = new1-new0; + + return kTRUE; +} + +/* +Bool_t Overlay(const TObject &obj) +{ + const TTL &ttl = static_cast(obj); + + const Double_t new0 = TMath::Max(fStart, ttl.fStart); + const Double_t new1 = TMath::Min(fStart+fLength, ttl.fStart+ttl.fLength); + + fStart = new0; + fLength = new1-new0; + + return IsValid(); +}*/ + + +void MDigitalSignal::Print(Option_t *o) const +{ + gLog << all << Form("%.2f,%.2f ", fStart, fStart+fLength) << endl; +} Index: trunk/MagicSoft/Mars/melectronics/MDigitalSignal.h =================================================================== --- trunk/MagicSoft/Mars/melectronics/MDigitalSignal.h (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/MDigitalSignal.h (revision 9242) @@ -0,0 +1,60 @@ +#ifndef MARS_MDigitalSignal +#define MARS_MDigitalSignal + +#ifndef ROOT_TObject +#include +#endif + +//#include + +class MDigitalSignal : public TObject +{ + Double_t fStart; + Double_t fLength; + +public: + MDigitalSignal(Double_t start=0, Double_t len=1e12/*FLT_MAX*/) : fStart(start), fLength(len) + { + } + + MDigitalSignal(const MDigitalSignal &ttl) : TObject(ttl), fStart(ttl.fStart), fLength(ttl.fLength) + { + } + + MDigitalSignal(const MDigitalSignal &ttl1, const MDigitalSignal &ttl2); + + // Getter + Bool_t IsValid() const { return fLength>0; } + + Double_t GetStart() const { return fStart; } + Double_t GetEnd() const { return fStart+fLength; } + Double_t GetLength() const { return fLength; } + + // TObject + Bool_t IsSortable() const { return kTRUE; } + Int_t Compare(const TObject *obj) const; + void Print(Option_t *o="") const; + + // MDigitalSignal + Bool_t Overlap(const TObject &obj) const; + Bool_t Combine(const TObject &obj); + + MDigitalSignal GetCoincidence(const MDigitalSignal &ttl) const + { + return MDigitalSignal(*this, ttl); + } + + const MDigitalSignal operator&&(const MDigitalSignal &ttl) const + { + return GetCoincidence(ttl); + } + const MDigitalSignal operator||(const MDigitalSignal &ttl) const + { + MDigitalSignal cpy(ttl); + return !cpy.Combine(*this) ? MDigitalSignal(0, -1) : cpy; + } + + ClassDef(MDigitalSignal, 1) // A digital signal with a start time and a length +}; + +#endif Index: trunk/MagicSoft/Mars/melectronics/Makefile =================================================================== --- trunk/MagicSoft/Mars/melectronics/Makefile (revision 9242) +++ trunk/MagicSoft/Mars/melectronics/Makefile (revision 9242) @@ -0,0 +1,34 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = Electronics + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mextralgo + +SRCFILES = MAvalanchePhotoDiode.cc \ + MAnalogSignal.cc \ + MAnalogChannels.cc \ + MDigitalSignal.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak