source: trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc@ 9521

Last change on this file since 9521 was 9518, checked in by tbretz, 15 years ago
*** empty log message ***
File size: 8.6 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of CheObs, the Modular Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appears in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: CheObs Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MSimAPD
28//
29// This tasks simulates the individual APDs. Before starting the APD is
30// initialized randomly according to the photon rate hitting the APD. Such
31// it is assumed that the initial condition of the APD is similar to the real
32// one. In this context it is assumed that the events are independent, so
33// that the APD is always in the same condition.
34//
35// For every photon and event the behaviour of the APD is simulated. The
36// output is set as weight to the MPhotonData containers.
37//
38// Remark:
39// - The photons MUST be sorted increasing in time.
40// - The photon rate used to initialize the APD must match the one used
41// to "fill" the random photons. (FIXME: This should be stored somewhere)
42//
43// Input Containers:
44// fNameGeomCam [MGeomCam]
45// MPhotonEvent
46// MPhotonStatistics
47//
48// Output Containers:
49// MPhotonEvent
50//
51//////////////////////////////////////////////////////////////////////////////
52#include "MSimAPD.h"
53
54#include <TH2.h>
55#include <TRandom.h>
56
57#include "MLog.h"
58#include "MLogManip.h"
59
60#include "MMath.h"
61#include "MParList.h"
62
63#include "MGeomCam.h"
64
65#include "MPhotonEvent.h"
66#include "MPhotonData.h"
67
68#include "MPedestalCam.h"
69#include "MPedestalPix.h"
70
71#include "MAvalanchePhotoDiode.h"
72
73ClassImp(MSimAPD);
74
75using namespace std;
76
77// --------------------------------------------------------------------------
78//
79// Default Constructor.
80//
81MSimAPD::MSimAPD(const char* name, const char *title)
82: fGeom(0), fEvt(0), fStat(0), fType(1)
83{
84 fName = name ? name : "MSimAPD";
85 fTitle = title ? title : " Task to simulate the detection behaviour of APDs";
86}
87
88// --------------------------------------------------------------------------
89//
90// Get the necessary parameter containers
91//
92Int_t MSimAPD::PreProcess(MParList *pList)
93{
94 if (fNameGeomCam.IsNull())
95 {
96 *fLog << inf << "No geometry container... skipping." << endl;
97 return kSKIP;
98 }
99
100 fGeom = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam");
101 if (!fGeom)
102 {
103 *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;
104
105 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
106 if (!fGeom)
107 {
108 *fLog << err << "MGeomCam not found... aborting." << endl;
109 return kFALSE;
110 }
111 }
112
113 fStat = (MPhotonStatistics*)pList->FindObject("MPhotonStatistics");
114 if (!fStat)
115 {
116 *fLog << err << "MPhotonStatistics not found... aborting." << endl;
117 return kFALSE;
118 }
119
120 fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
121 if (!fEvt)
122 {
123 *fLog << err << "MPhotonEvent not found... aborting." << endl;
124 return kFALSE;
125 }
126
127 fRates = (MPedestalCam*)pList->FindObject("AccidentalPhotonRates", "MPedestalCam");
128 if (!fRates)
129 {
130 *fLog << inf;
131 *fLog << "AccidentalPhotonRates [MPedestalCam] not found..." << endl;
132 *fLog << " using " << fFreq << " as default for all G-APDs." << endl;
133 }
134
135 return kTRUE;
136}
137
138// --------------------------------------------------------------------------
139//
140// Initialize as many APDs as we have pixels in the fGeomCam
141//
142Bool_t MSimAPD::ReInit(MParList *plist)
143{
144 if (UInt_t(fAPDs.GetEntriesFast())==fGeom->GetNumPixels())
145 return kTRUE;
146
147 fAPDs.Delete();
148
149 // FIXME:
150 // * initialize an empty APD and read the APD setup from a file to
151 // allow different APDs.
152 // * Make the arguments a data member of MSimAPD
153
154 Int_t ncells = 0;
155 Float_t crosstalk = 0;
156 Float_t deadtime = 0;
157 Float_t recovery = 0;
158
159 switch (fType)
160 {
161 case 1:
162 ncells = 30;
163 crosstalk = 0.2;
164 deadtime = 3;
165 recovery = 8.75*4;
166 break;
167
168 case 2:
169 ncells = 60;
170 crosstalk = 0.2;
171 deadtime = 3;
172 recovery = 8.75;
173 break;
174
175 case 3:
176 ncells = 60;
177 crosstalk = 0.15;
178 deadtime = 3;
179 recovery = 8.75;
180 break;
181
182 default:
183 *fLog << err << "ERROR - APD type " << fType << " undefined." << endl;
184 return kFALSE;
185 }
186
187 for (UInt_t i=0; i<fGeom->GetNumPixels(); i++)
188 fAPDs.Add(new APD(ncells, crosstalk, deadtime, recovery));
189
190 return kTRUE;
191}
192
193// --------------------------------------------------------------------------
194//
195// Process all photons through the corresponding APD and set the output
196// (weight) accordingly.
197//
198Int_t MSimAPD::Process()
199{
200 //const Double_t rate = 40e9;
201 // FIXME: Where do we get this number from??
202 // const Double_t rate = 0.04;
203
204 // Make all APDs look neutral for the first hit by a photon according to the
205 // average hit rate
206 const UInt_t npix = fAPDs.GetEntriesFast();
207
208 // Check if we can safely proceed (this can fail if we either haven't been
209 // ReInit'ed or the max index in MPhotonStatistics is wrong)
210 if ((Int_t)npix<fStat->GetMaxIndex())
211 {
212 *fLog << err << "ERROR - MSimAPD::Process: Only " << npix << " APDs initialized. At least " << fStat->GetMaxIndex() << " needed... abort." << endl;
213 return kERROR;
214 }
215/*
216 for (UInt_t idx=0; idx<npix; idx++)
217 {
218 const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
219 static_cast<APD*>(fAPDs.UncheckedAt(idx))->FillRandom(freq, fStat->GetTimeFirst());
220 }
221*/
222
223 // This tries to initialize dead and relaxing cells properly. If
224 // the APD has not been initialized before the chip is randomsly
225 // filled, otherwise a time window of the default relaxing time
226 // is simulated, so that the previous influence is less than a permille.
227 for (UInt_t idx=0; idx<npix; idx++)
228 {
229 const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
230 static_cast<APD*>(fAPDs.UncheckedAt(idx))->Init(freq);
231 }
232
233 // Get number of photons
234 const Int_t num = fEvt->GetNumPhotons();
235
236 // Loop over all photons
237 for (Int_t i=0; i<num; i++)
238 {
239 // Get i-th photon
240 MPhotonData &ph = (*fEvt)[i];
241
242 // Get arrival time of photon and idx
243 const Double_t t = ph.GetTime()-fStat->GetTimeFirst();
244 const Int_t idx = ph.GetTag();
245 if (idx<0)
246 {
247 *fLog << err << "ERROR - MSimAPD: Invalid index -1." << endl;
248 return kERROR;
249 }
250
251 if (ph.GetWeight()!=1)
252 {
253 *fLog << err << "ERROR - MSimAPD: Weight of " << i << "-th photon not 1, but " << ph.GetWeight() << endl;
254 ph.Print();
255 return kERROR;
256 }
257 // Simulate hitting the APD (the signal height in effective
258 // "number of photons" is returned)
259 const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCellRelative(t);
260
261 // FIXME: Make a proper simulation of the excess noise!!!
262 //const Double_t signal = gRandom->Gaus(hits, 0.2*TMath::Sqrt(hits));
263
264 // Set the weight to the input
265 ph.SetWeight(hits);
266 }
267
268 // Now we have to shift the evolved time of all APDs to the end of our
269 // simulated time.
270 for (UInt_t idx=0; idx<npix; idx++)
271 static_cast<APD*>(fAPDs.UncheckedAt(idx))->IncreaseTime(fStat->GetTimeLast());
272
273 return kTRUE;
274}
275
276// --------------------------------------------------------------------------
277//
278// NameGeomCam
279// Type: 1
280//
281Int_t MSimAPD::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
282{
283 Bool_t rc = kFALSE;
284 if (IsEnvDefined(env, prefix, "NameGeomCam", print))
285 {
286 rc = kTRUE;
287 fNameGeomCam = GetEnvValue(env, prefix, "NameGeomCam", fNameGeomCam);
288 }
289
290 if (IsEnvDefined(env, prefix, "Type", print))
291 {
292 rc = kTRUE;
293 fType = GetEnvValue(env, prefix, "Type", fType);
294 }
295
296 return rc;
297}
Note: See TracBrowser for help on using the repository browser.