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 | // MSimTrigger
|
---|
28 | //
|
---|
29 | // This task takes the pure analog channels and simulates a trigger
|
---|
30 | // electronics.
|
---|
31 | //
|
---|
32 | // In a first step several channels can be summed together by a look-up table
|
---|
33 | // fRouteAC.
|
---|
34 | //
|
---|
35 | // In a second step from these analog channels the output of a discriminator
|
---|
36 | // is calculated using a threshold and optional a fixed digital signal length.
|
---|
37 | //
|
---|
38 | // The signal length of the digital signal emitted by the discriminator
|
---|
39 | // can either be bound to the time the signal is above the threshold
|
---|
40 | // defined by fDiscriminatorThreshold if fDigitalSignalLength<0 or set to a
|
---|
41 | // fixed length (fDigitalSignalLength>0).
|
---|
42 | //
|
---|
43 | // With a second look-up table fCoincidenceMap the analog channels are
|
---|
44 | // checked for coincidences. The coincidence must at least be of the length
|
---|
45 | // defined by fCoincidenceTime. The earliest coincide is then stored as
|
---|
46 | // trigger position.
|
---|
47 | //
|
---|
48 | // If a minimum multiplicity m is given, m signals above threshold
|
---|
49 | // in the coincidence patterns are enough to emit a trigger signal.
|
---|
50 | //
|
---|
51 | //
|
---|
52 | // For MAGIC1:
|
---|
53 | // - fDigitalSignalLength between 6ns and 12ns
|
---|
54 | // - fCoincidenceTime between 0.25ns to 1ns
|
---|
55 | //
|
---|
56 | //
|
---|
57 | // Input Containers:
|
---|
58 | // IntendedPulsePos [MParameterD]
|
---|
59 | // MAnalogChannels
|
---|
60 | // MRawRunHeader
|
---|
61 | //
|
---|
62 | // Output Containers:
|
---|
63 | // TriggerPos [MParameterD]
|
---|
64 | // MRawEvtHeader
|
---|
65 | //
|
---|
66 | //////////////////////////////////////////////////////////////////////////////
|
---|
67 | #include "MSimTrigger.h"
|
---|
68 |
|
---|
69 | #include "MLog.h"
|
---|
70 | #include "MLogManip.h"
|
---|
71 |
|
---|
72 | #include "MParList.h"
|
---|
73 | #include "MParameters.h"
|
---|
74 |
|
---|
75 | #include "MLut.h"
|
---|
76 | #include "MArrayI.h"
|
---|
77 |
|
---|
78 | #include "MRawEvtHeader.h"
|
---|
79 | #include "MRawRunHeader.h"
|
---|
80 |
|
---|
81 | #include "MAnalogSignal.h"
|
---|
82 | #include "MAnalogChannels.h"
|
---|
83 | #include "MDigitalSignal.h"
|
---|
84 |
|
---|
85 | #include "MTriggerPattern.h"
|
---|
86 |
|
---|
87 | #include "MPedestalCam.h"
|
---|
88 | #include "MPedestalPix.h"
|
---|
89 |
|
---|
90 | ClassImp(MSimTrigger);
|
---|
91 |
|
---|
92 | using namespace std;
|
---|
93 |
|
---|
94 | // --------------------------------------------------------------------------
|
---|
95 | //
|
---|
96 | // Default Constructor.
|
---|
97 | //
|
---|
98 | MSimTrigger::MSimTrigger(const char *name, const char *title)
|
---|
99 | : fCamera(0),
|
---|
100 | fPulsePos(0),
|
---|
101 | fTrigger(0),
|
---|
102 | fRunHeader(0),
|
---|
103 | fEvtHeader(0),
|
---|
104 | fElectronicNoise(0),
|
---|
105 | fGain(0),
|
---|
106 | fDiscriminatorThreshold(-1),
|
---|
107 | fDigitalSignalLength(8),
|
---|
108 | fCoincidenceTime(0.5),
|
---|
109 | fShiftBaseline(kTRUE),
|
---|
110 | fUngainSignal(kTRUE),
|
---|
111 | fSimulateElectronics(kTRUE),
|
---|
112 | fMinMultiplicity(-1),
|
---|
113 | fCableDelay(21),
|
---|
114 | fCableDamping(0.) // default Damping Set to zero, so users, who do not set
|
---|
115 | // the CableDamoing in the ceres.rc do not see a difference.
|
---|
116 |
|
---|
117 | {
|
---|
118 | fName = name ? name : "MSimTrigger";
|
---|
119 | fTitle = title ? title : "Task to simulate trigger electronics";
|
---|
120 | }
|
---|
121 |
|
---|
122 | // --------------------------------------------------------------------------
|
---|
123 | //
|
---|
124 | // Take two TObjArrays with a collection of digital signals.
|
---|
125 | // Every signal from one array is compared with any from the other array.
|
---|
126 | // For all signals which overlap and which have an overlap time >gate
|
---|
127 | // a new digital signal is created storing start time and length of overlap.
|
---|
128 | // They are collected in a newly allocated TObjArray. A pointer to this array
|
---|
129 | // is returned.
|
---|
130 | //
|
---|
131 | // The user gains owenership of the object, i.e., the user is responsible of
|
---|
132 | // deleting the memory.
|
---|
133 | //
|
---|
134 | TObjArray *MSimTrigger::CalcCoincidence(const TObjArray &arr1, const TObjArray &arr2/*, Float_t gate*/) const
|
---|
135 | {
|
---|
136 | TObjArray *res = new TObjArray;
|
---|
137 |
|
---|
138 | if (arr1.GetEntriesFast()==0 || arr2.GetEntriesFast()==0)
|
---|
139 | return res;
|
---|
140 |
|
---|
141 | TIter Next1(&arr1);
|
---|
142 | MDigitalSignal *ttl1 = 0;
|
---|
143 | while ((ttl1=static_cast<MDigitalSignal*>(Next1())))
|
---|
144 | {
|
---|
145 | TIter Next2(&arr2);
|
---|
146 | MDigitalSignal *ttl2 = 0;
|
---|
147 | while ((ttl2=static_cast<MDigitalSignal*>(Next2())))
|
---|
148 | {
|
---|
149 | MDigitalSignal *ttl = new MDigitalSignal(*ttl1, *ttl2);
|
---|
150 | /*
|
---|
151 | if (ttl->GetLength()<=gate)
|
---|
152 | {
|
---|
153 | delete ttl;
|
---|
154 | continue;
|
---|
155 | }
|
---|
156 | */
|
---|
157 | res->Add(ttl);
|
---|
158 | }
|
---|
159 | }
|
---|
160 |
|
---|
161 | res->SetOwner();
|
---|
162 |
|
---|
163 | return res;
|
---|
164 | }
|
---|
165 |
|
---|
166 | class Edge : public TObject
|
---|
167 | {
|
---|
168 | private:
|
---|
169 | Double_t fEdge;
|
---|
170 | Int_t fRising;
|
---|
171 |
|
---|
172 | public:
|
---|
173 | Edge(Double_t t, Int_t rising) : fEdge(t), fRising(rising) { }
|
---|
174 | Bool_t IsSortable() const { return kTRUE; }
|
---|
175 | Int_t Compare(const TObject *o) const { const Edge *e = static_cast<const Edge*>(o); if (e->fEdge<fEdge) return 1; if (e->fEdge>fEdge) return -1; return 0; }
|
---|
176 |
|
---|
177 | Int_t IsRising() const { return fRising; }
|
---|
178 | Double_t GetEdge() const { return fEdge; }
|
---|
179 | };
|
---|
180 |
|
---|
181 | // --------------------------------------------------------------------------
|
---|
182 | //
|
---|
183 | // Calculate a multiplicity trigger on the given array(s). The idx-array
|
---|
184 | // conatins all channels which should be checked for coincidences
|
---|
185 | // and the ttls array conatins the arrays with the digital signals.
|
---|
186 | //
|
---|
187 | // For the windows in which more or euqal than threshold channels have
|
---|
188 | // a high signal a new MDigitalSignal is created. newly allocated
|
---|
189 | // array with a collection of these trigger signals is returned.
|
---|
190 | //
|
---|
191 | TObjArray *MSimTrigger::CalcMinMultiplicity(const MArrayI &idx, const TObjArray &ttls, Int_t threshold) const
|
---|
192 | {
|
---|
193 | // Create a new array for the rising and falling edges of the signals
|
---|
194 | TObjArray times;
|
---|
195 | times.SetOwner();
|
---|
196 |
|
---|
197 | // Fill the array with edges from all digital signals of all our channels
|
---|
198 | for (UInt_t k=0; k<idx.GetSize(); k++)
|
---|
199 | {
|
---|
200 | TObjArray *arr = static_cast<TObjArray*>(ttls[idx[k]]);
|
---|
201 |
|
---|
202 | TIter Next(arr);
|
---|
203 | MDigitalSignal *ttl = 0;
|
---|
204 | while ((ttl=static_cast<MDigitalSignal*>(Next())))
|
---|
205 | {
|
---|
206 | times.Add(new Edge(ttl->GetStart(), 1));
|
---|
207 | times.Add(new Edge(ttl->GetEnd(), -1));
|
---|
208 | }
|
---|
209 | }
|
---|
210 |
|
---|
211 | // Sort them in time
|
---|
212 | times.Sort();
|
---|
213 |
|
---|
214 | // Start with no channel active
|
---|
215 | Int_t lvl = 0;
|
---|
216 |
|
---|
217 | TObjArray *res = new TObjArray;
|
---|
218 | res->SetOwner();
|
---|
219 |
|
---|
220 | // First remove all edges which do not change the status
|
---|
221 | // "below threshold" or "above threshold"
|
---|
222 | for (int i=0; i<times.GetEntriesFast(); i++)
|
---|
223 | {
|
---|
224 | // Get i-th edge
|
---|
225 | const Edge &e = *static_cast<Edge*>(times.UncheckedAt(i));
|
---|
226 |
|
---|
227 | // Claculate what the number of active channels after the edge is
|
---|
228 | const Int_t lvl1 = lvl + e.IsRising();
|
---|
229 |
|
---|
230 | // Remove edge if number of active channels before and after the
|
---|
231 | // edge lower is lower than the threshold or higher than
|
---|
232 | // the threshold
|
---|
233 | if (lvl+1<threshold || lvl-1>=threshold)
|
---|
234 | delete times.RemoveAt(i);
|
---|
235 |
|
---|
236 | // keep the (now) "previous" level
|
---|
237 | lvl = lvl1<0 ? 0 : lvl1;
|
---|
238 | }
|
---|
239 |
|
---|
240 | // Remove the empty slots from the array
|
---|
241 | times.Compress();
|
---|
242 |
|
---|
243 | //
|
---|
244 | for (int i=0; i<times.GetEntriesFast()-1; i++)
|
---|
245 | {
|
---|
246 | // get the current edge
|
---|
247 | const Edge &e0 = *static_cast<Edge*>(times.UncheckedAt(i));
|
---|
248 |
|
---|
249 | // go ahead if this is a falling edge
|
---|
250 | if (e0.IsRising()!=1)
|
---|
251 | continue;
|
---|
252 |
|
---|
253 | // get the following edge (must be a falling edge now)
|
---|
254 | const Edge &e1 = *static_cast<Edge*>(times.UncheckedAt(i+1));
|
---|
255 |
|
---|
256 | // calculate the length of the digital signal
|
---|
257 | const Double_t len = e1.GetEdge()-e0.GetEdge();
|
---|
258 |
|
---|
259 | // Create a digital trigger signal
|
---|
260 | MDigitalSignal *ds = new MDigitalSignal(e0.GetEdge(), len);
|
---|
261 | //ds->SetIndex(lvl);
|
---|
262 | res->Add(ds);
|
---|
263 | }
|
---|
264 |
|
---|
265 | return res;
|
---|
266 | }
|
---|
267 |
|
---|
268 | // --------------------------------------------------------------------------
|
---|
269 | //
|
---|
270 | // Check for the necessary parameter containers. Read the luts.
|
---|
271 | //
|
---|
272 | Int_t MSimTrigger::PreProcess(MParList *pList)
|
---|
273 | {
|
---|
274 | fTrigger = (MParameterD*)pList->FindCreateObj("MParameterD", "TriggerPos");
|
---|
275 | if (!fTrigger)
|
---|
276 | return kFALSE;
|
---|
277 |
|
---|
278 | fPulsePos = (MParameterD*)pList->FindObject("IntendedPulsePos", "MParameterD");
|
---|
279 | if (!fPulsePos)
|
---|
280 | {
|
---|
281 | *fLog << err << "IntendedPulsePos [MParameterD] not found... aborting." << endl;
|
---|
282 | return kFALSE;
|
---|
283 | }
|
---|
284 |
|
---|
285 | fCamera = (MAnalogChannels*)pList->FindObject("MAnalogChannels");
|
---|
286 | if (!fCamera)
|
---|
287 | {
|
---|
288 | *fLog << err << "MAnalogChannels not found... aborting." << endl;
|
---|
289 | return kFALSE;
|
---|
290 | }
|
---|
291 |
|
---|
292 | fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
|
---|
293 | if (!fRunHeader)
|
---|
294 | {
|
---|
295 | *fLog << err << "MRawRunHeader not found... aborting." << endl;
|
---|
296 | return kFALSE;
|
---|
297 | }
|
---|
298 |
|
---|
299 | fEvtHeader = (MRawEvtHeader*)pList->FindCreateObj("MRawEvtHeader");
|
---|
300 | if (!fEvtHeader)
|
---|
301 | return kFALSE;
|
---|
302 |
|
---|
303 | if (!fSimulateElectronics)
|
---|
304 | {
|
---|
305 | *fLog << inf << "Simulation of electronics switched off... first photon will trigger." << endl;
|
---|
306 | return kTRUE;
|
---|
307 | }
|
---|
308 |
|
---|
309 | fElectronicNoise = 0;
|
---|
310 | if (fShiftBaseline)
|
---|
311 | {
|
---|
312 | fElectronicNoise = (MPedestalCam*)pList->FindObject("ElectronicNoise", "MPedestalCam");
|
---|
313 | if (!fElectronicNoise)
|
---|
314 | {
|
---|
315 | *fLog << err << "ElectronicNoise [MPedestalCam] not found... aborting." << endl;
|
---|
316 | return kFALSE;
|
---|
317 | }
|
---|
318 | *fLog << inf << "Baseline will be shifted back to 0 for discriminator." << endl;
|
---|
319 | }
|
---|
320 |
|
---|
321 | fGain = 0;
|
---|
322 | if (fUngainSignal)
|
---|
323 | {
|
---|
324 | fGain = (MPedestalCam*)pList->FindObject("Gain", "MPedestalCam");
|
---|
325 | if (!fGain)
|
---|
326 | {
|
---|
327 | *fLog << err << "Gain [MPedestalCam] not found... aborting." << endl;
|
---|
328 | return kFALSE;
|
---|
329 | }
|
---|
330 | *fLog << inf << "Discriminator will be multiplied by applied gain." << endl;
|
---|
331 | }
|
---|
332 |
|
---|
333 | fRouteAC.Delete();
|
---|
334 | if (!fNameRouteAC.IsNull() && fRouteAC.ReadFile(fNameRouteAC)<0)
|
---|
335 | return kFALSE;
|
---|
336 |
|
---|
337 | fCoincidenceMap.Delete();
|
---|
338 | if (!fNameCoincidenceMap.IsNull() && fCoincidenceMap.ReadFile(fNameCoincidenceMap)<0)
|
---|
339 | return kFALSE;
|
---|
340 |
|
---|
341 | // ---------------- Consistency checks ----------------------
|
---|
342 |
|
---|
343 | if (!fRouteAC.IsEmpty() && !fCoincidenceMap.IsEmpty() &&
|
---|
344 | fCoincidenceMap.GetMaxIndex()>fRouteAC.GetNumRows()-1)
|
---|
345 | {
|
---|
346 | *fLog << err;
|
---|
347 | *fLog << "ERROR - AC routing produces " << fRouteAC.GetNumRows() << " analog channels," << endl;
|
---|
348 | *fLog << " but the coincidence map expects at least " << fCoincidenceMap.GetMaxIndex()+1 << " channels." << endl;
|
---|
349 | return kERROR;
|
---|
350 | }
|
---|
351 |
|
---|
352 | // if (fDiscriminatorThreshold<=0)
|
---|
353 | // {
|
---|
354 | // *fLog << err << "ERROR - Discriminator threshold " << fDiscriminatorThreshold << " invalid." << endl;
|
---|
355 | // return kFALSE;
|
---|
356 | // }
|
---|
357 |
|
---|
358 | if (fElectronicNoise && !fRouteAC.IsEmpty() && !fRouteAC.IsDefaultCol())
|
---|
359 | {
|
---|
360 | // FIXME: Apply to analog channels when summing
|
---|
361 | *fLog << warn << "WARNING - A baseline shift doesn't make sense for sum-channels... reset." << endl;
|
---|
362 | fElectronicNoise = 0;
|
---|
363 | }
|
---|
364 |
|
---|
365 | if (fGain && !fRouteAC.IsEmpty() && !fRouteAC.IsDefaultCol())
|
---|
366 | {
|
---|
367 | // FIXME: Apply to analog channels when summing
|
---|
368 | *fLog << warn << "WARNING - Ungain doesn't make sense for sum-channels... reset." << endl;
|
---|
369 | fGain = 0;
|
---|
370 | }
|
---|
371 |
|
---|
372 |
|
---|
373 | // ---------------- Information output ----------------------
|
---|
374 |
|
---|
375 | *fLog << inf;
|
---|
376 |
|
---|
377 | if (fRouteAC.IsEmpty())
|
---|
378 | *fLog << "Re-routing/summing of analog channels before discriminator switched off." << endl;
|
---|
379 | else
|
---|
380 | *fLog << "Using " << fNameRouteAC << " for re-routing/summing of analog channels before discriminator." << endl;
|
---|
381 |
|
---|
382 | if (fCoincidenceMap.IsEmpty() && fMinMultiplicity<=0)
|
---|
383 | *fLog << "No coincidences of digital channels will be checked. Signal-above-threshold trigger applied." << endl;
|
---|
384 | else
|
---|
385 | {
|
---|
386 | *fLog << "Using ";
|
---|
387 | if (fCoincidenceMap.IsEmpty())
|
---|
388 | *fLog << "the whole camera";
|
---|
389 | else
|
---|
390 | *fLog << "patterns from " << fNameCoincidenceMap;
|
---|
391 | *fLog << " to check for ";
|
---|
392 | if (fMinMultiplicity>0)
|
---|
393 | *fLog << fMinMultiplicity << " multiplicity." << endl;
|
---|
394 | else
|
---|
395 | *fLog << "coincidences of the digital channels." << endl;
|
---|
396 | }
|
---|
397 |
|
---|
398 | *fLog << "Using discriminator threshold of " << fDiscriminatorThreshold << endl;
|
---|
399 |
|
---|
400 | *fLog << "Using fCableDelay " << fCableDelay << endl;
|
---|
401 | *fLog << "Using fCableDamping " << fCableDamping << endl;
|
---|
402 |
|
---|
403 | ////////////////////////////////
|
---|
404 | // open some output files for debugging
|
---|
405 | // patch_file.open("/home/fact_opr/patch_file.csv", ios_base::out);
|
---|
406 | // clipped_file.open("/home/fact_opr/clipped_file.csv", ios_base::out);
|
---|
407 | // digital_file.open("/home/fact_opr/digital_file.csv", ios_base::out);
|
---|
408 | // ratescan_file.open("/home/fact_opr/ratescan_file.csv", ios_base::out);
|
---|
409 |
|
---|
410 |
|
---|
411 |
|
---|
412 | return kTRUE;
|
---|
413 | }
|
---|
414 |
|
---|
415 | /*
|
---|
416 | class MDigitalChannel : public TObjArray
|
---|
417 | {
|
---|
418 | private:
|
---|
419 | TObjArray fArray;
|
---|
420 |
|
---|
421 | public:
|
---|
422 | MDigitalSignal *GetSignal(UInt_t i) { return static_cast<MDigitalSignal*>(fArray[i]); }
|
---|
423 |
|
---|
424 | };
|
---|
425 | */
|
---|
426 |
|
---|
427 | #include "MCamEvent.h"
|
---|
428 | class MTriggerSignal : public MParContainer, public MCamEvent
|
---|
429 | {
|
---|
430 | private:
|
---|
431 | TObjArray fSignals;
|
---|
432 |
|
---|
433 | public:
|
---|
434 | MTriggerSignal() { fSignals.SetOwner(); }
|
---|
435 |
|
---|
436 | void Add(MDigitalSignal *signal) { fSignals.Add(signal); }
|
---|
437 |
|
---|
438 | MDigitalSignal *GetSignal(UInt_t i) { return static_cast<MDigitalSignal*>(fSignals[i]); }
|
---|
439 |
|
---|
440 | void Sort() { fSignals.Sort(); }
|
---|
441 |
|
---|
442 | Int_t GetNumSignals() const { return fSignals.GetEntriesFast(); }
|
---|
443 |
|
---|
444 | Float_t GetFirstTrigger() const
|
---|
445 | {
|
---|
446 | MDigitalSignal *sig = static_cast<MDigitalSignal*>(fSignals[0]);
|
---|
447 | return sig ? sig->GetStart() : -1;
|
---|
448 | }
|
---|
449 | Int_t GetFirstIndex() const
|
---|
450 | {
|
---|
451 | MDigitalSignal *sig = static_cast<MDigitalSignal*>(fSignals[0]);
|
---|
452 | return sig ? sig->GetIndex() : -1;
|
---|
453 | }
|
---|
454 | Bool_t GetPixelContent(Double_t&, Int_t, const MGeomCam&, Int_t) const
|
---|
455 | {
|
---|
456 | switch (1)
|
---|
457 | {
|
---|
458 | case 1: // yes/no
|
---|
459 | case 2: // first time
|
---|
460 | case 3: // length
|
---|
461 | case 4: // n
|
---|
462 | break;
|
---|
463 | }
|
---|
464 |
|
---|
465 | return kTRUE;
|
---|
466 | }
|
---|
467 | void DrawPixelContent(Int_t) const
|
---|
468 | {
|
---|
469 | }
|
---|
470 | };
|
---|
471 |
|
---|
472 |
|
---|
473 | void MSimTrigger::SetTrigger(Double_t pos, Int_t idx)
|
---|
474 | {
|
---|
475 | // FIXME: Jitter! (Own class?)
|
---|
476 | fTrigger->SetVal(pos);
|
---|
477 | fTrigger->SetReadyToSave();
|
---|
478 |
|
---|
479 | // Flag this event as triggered by the lvl1 trigger (FIXME?)
|
---|
480 | fEvtHeader->SetTriggerPattern(pos<0 ? 0 : 4);
|
---|
481 | fEvtHeader->SetNumTrigLvl1((UInt_t)idx);
|
---|
482 | fEvtHeader->SetReadyToSave();
|
---|
483 | }
|
---|
484 |
|
---|
485 | // --------------------------------------------------------------------------
|
---|
486 | //
|
---|
487 | Int_t MSimTrigger::Process()
|
---|
488 | {
|
---|
489 | // Invalidate trigger
|
---|
490 | //fTrigger->SetVal(-1);
|
---|
491 | // Calculate the minimum and maximum time for a valid trigger
|
---|
492 | const Double_t freq = fRunHeader->GetFreqSampling()/1000.;
|
---|
493 | const Float_t nsamp = fRunHeader->GetNumSamplesHiGain();
|
---|
494 | const Float_t pulspos = fPulsePos->GetVal()/freq;
|
---|
495 |
|
---|
496 | // Valid range in units of bins
|
---|
497 | const Float_t min = fCamera->GetValidRangeMin()+pulspos;
|
---|
498 | const Float_t max = fCamera->GetValidRangeMax()-(nsamp-pulspos);
|
---|
499 |
|
---|
500 | if (!fSimulateElectronics)
|
---|
501 | {
|
---|
502 | SetTrigger(min, -1);
|
---|
503 | return kTRUE;
|
---|
504 | }
|
---|
505 |
|
---|
506 | // ================== Simulate channel bundling ====================
|
---|
507 |
|
---|
508 | // FIXME: Before we can bundle the channels we have to make a copy
|
---|
509 | // and simulate clipping
|
---|
510 |
|
---|
511 | // Check if routing should be done
|
---|
512 | const Bool_t empty = fRouteAC.IsEmpty();
|
---|
513 |
|
---|
514 | // If no channels are summed the number of patches stays the same
|
---|
515 | const UInt_t npatch = empty ? fCamera->GetNumChannels() : fRouteAC.GetEntriesFast();
|
---|
516 |
|
---|
517 | // Use the given analog channels as default out. If channels are
|
---|
518 | // summed overwrite with a newly allocated set of analog channels
|
---|
519 | MAnalogChannels *patches = fCamera;
|
---|
520 | MAnalogChannels *raw_patches = fCamera;
|
---|
521 | if (!empty)
|
---|
522 | {
|
---|
523 | // FIXME: Can we add gain and offset here into a new container?
|
---|
524 |
|
---|
525 | patches = new MAnalogChannels(npatch, fCamera->GetNumSamples());
|
---|
526 | raw_patches = new MAnalogChannels(npatch, fCamera->GetNumSamples());
|
---|
527 | for (UInt_t patch_id=0; patch_id<npatch; patch_id++)
|
---|
528 | {
|
---|
529 | const MArrayI &row = fRouteAC.GetRow(patch_id);
|
---|
530 | for (UInt_t pxl_in_patch=0; pxl_in_patch<row.GetSize(); pxl_in_patch++)
|
---|
531 | {
|
---|
532 | const UInt_t pixel_id = row[pxl_in_patch];
|
---|
533 |
|
---|
534 | // FIXME: Shrinking the mapping table earlier (e.g.
|
---|
535 | // ReInit) would avoid a lot of if's
|
---|
536 | if (pixel_id<fCamera->GetNumChannels())
|
---|
537 | {
|
---|
538 | (*raw_patches)[patch_id].AddSignal((*fCamera)[pixel_id]);
|
---|
539 | (*patches)[patch_id].AddSignal((*fCamera)[pixel_id]);
|
---|
540 | (*patches)[patch_id].AddSignal((*fCamera)[pixel_id], fCableDelay, fCableDamping);
|
---|
541 | }
|
---|
542 | }
|
---|
543 | }
|
---|
544 | }
|
---|
545 |
|
---|
546 | // DN: 20140219 Ratescan:
|
---|
547 | //
|
---|
548 | //
|
---|
549 | // for (UInt_t patch_id=0; patch_id<npatch; patch_id++)
|
---|
550 | // {
|
---|
551 | // MAnalogSignal current_patch = (*raw_patches)[patch_id];
|
---|
552 | // float max = current_patch[0];
|
---|
553 | // for (UInt_t i=1; i<current_patch.GetSize(); i++)
|
---|
554 | // {
|
---|
555 | // if (current_patch[i] > max)
|
---|
556 | // {
|
---|
557 | // max = current_patch[i];
|
---|
558 | // }
|
---|
559 | // }
|
---|
560 | // ratescan_file << max << " ";
|
---|
561 | // }
|
---|
562 | // ratescan_file << endl;
|
---|
563 |
|
---|
564 | // // DN 20131108: DEBUGGING:
|
---|
565 | // for (UInt_t patch_id=0; patch_id<npatch; patch_id++)
|
---|
566 | // {
|
---|
567 | // for (UInt_t slice=0; slice<fCamera->GetNumSamples(); slice++)
|
---|
568 | // {
|
---|
569 | // patch_file << (*raw_patches)[patch_id][slice] << "\t";
|
---|
570 | // clipped_file << (*patches)[patch_id][slice] << "\t";
|
---|
571 | // }
|
---|
572 | // patch_file << endl;
|
---|
573 | // clipped_file << endl;
|
---|
574 | // }
|
---|
575 |
|
---|
576 |
|
---|
577 |
|
---|
578 | // FIXME: Write patches
|
---|
579 |
|
---|
580 | // ================== Simulate discriminators ====================
|
---|
581 |
|
---|
582 | TObjArray ttls(npatch);
|
---|
583 | ttls.SetOwner();
|
---|
584 |
|
---|
585 | for (UInt_t i=0; i<npatch; i++)
|
---|
586 | {
|
---|
587 | // FIXME: What if the gain was also allpied to the baseline?
|
---|
588 | const Double_t offset = fElectronicNoise ? (*fElectronicNoise)[i].GetPedestal() : 0;
|
---|
589 | const Double_t gain = fGain ? (*fGain)[i].GetPedestal() : 1;
|
---|
590 | ttls.AddAt(
|
---|
591 | (*patches)[i].Discriminate(
|
---|
592 | fDiscriminatorThreshold*gain+offset, // treshold
|
---|
593 | Double_t(fCableDelay), // start
|
---|
594 | Double_t(fCamera->GetNumSamples() - fCableDelay), // end
|
---|
595 | //fDigitalSignalLength // time-over-threshold, or fixed-length?
|
---|
596 | -1 // -1 = time-over-threshold
|
---|
597 | ),
|
---|
598 | i
|
---|
599 | );
|
---|
600 | }
|
---|
601 |
|
---|
602 | // FIXME: Write TTLs!
|
---|
603 |
|
---|
604 | // If analog channels had been newly allocated free memmory
|
---|
605 | if (patches!=fCamera)
|
---|
606 | delete patches;
|
---|
607 | if (raw_patches!=fCamera)
|
---|
608 | delete raw_patches;
|
---|
609 |
|
---|
610 | // =================== Simulate coincidences ======================
|
---|
611 |
|
---|
612 | // If the map is empty we create a one-pixel-coincidence map
|
---|
613 | // FIMXE: This could maybe be accelerated if the Clone can be
|
---|
614 | // omitted in the loop
|
---|
615 | if (fCoincidenceMap.IsEmpty())
|
---|
616 | {
|
---|
617 | if (fMinMultiplicity>0)
|
---|
618 | fCoincidenceMap.SetDefaultRow(npatch);
|
---|
619 | else
|
---|
620 | fCoincidenceMap.SetDefaultCol(npatch);
|
---|
621 | }
|
---|
622 |
|
---|
623 | // Create an array for the individual triggers
|
---|
624 | MTriggerSignal triggers;
|
---|
625 |
|
---|
626 | Int_t cnt = 0;
|
---|
627 | Int_t rmlo = 0;
|
---|
628 | Int_t rmhi = 0;
|
---|
629 |
|
---|
630 | // cout << "MSimTrigger::fMinMultiplicity = " << fMinMultiplicity << endl;
|
---|
631 | // cout << "MSimTrigger::fCoincidenceTime = " << fCoincidenceTime << endl;
|
---|
632 | // cout << "fCoincidenceMap.GetEntries() = " << fCoincidenceMap.GetEntries() << endl;
|
---|
633 | // cout << "MSimTrigger::fCableDelay = " << fCableDelay << endl;
|
---|
634 | // cout << "MSimTrigger::fCableDamping = " << fCableDamping << endl;
|
---|
635 | // cout << "min:" << min << endl;
|
---|
636 | // cout << "max:" << max << endl;
|
---|
637 |
|
---|
638 | for (int j=0; j<fCoincidenceMap.GetEntries(); j++)
|
---|
639 | {
|
---|
640 | const MArrayI &idx = fCoincidenceMap.GetRow(j);
|
---|
641 |
|
---|
642 | TObjArray *arr = 0;
|
---|
643 |
|
---|
644 | if (fMinMultiplicity>0)
|
---|
645 | {
|
---|
646 | arr = CalcMinMultiplicity(idx, ttls, fMinMultiplicity);
|
---|
647 | }
|
---|
648 | else
|
---|
649 | {
|
---|
650 | arr = CalcMinMultiplicity(idx, ttls, idx.GetSize());
|
---|
651 | /*
|
---|
652 | // Start with a copy of the first coincidence channel
|
---|
653 | arr = static_cast<TObjArray*>(ttls[idx[0]]->Clone());
|
---|
654 | arr->SetOwner();
|
---|
655 |
|
---|
656 | // compare to all other channels in this coincidence patch, one by one
|
---|
657 | for (UInt_t k=1; k<idx.GetSize() && arr->GetEntriesFast()>0; k++)
|
---|
658 | {
|
---|
659 | TObjArray *res = CalcCoincidence(*arr, *static_cast<TObjArray*>(ttls[idx[k]]));//, fCoincidenceTime);
|
---|
660 |
|
---|
661 | // Delete the original array and keep the new one
|
---|
662 | delete arr;
|
---|
663 | arr = res;
|
---|
664 | }*/
|
---|
665 | }
|
---|
666 |
|
---|
667 | // Count the number of totally emitted coincidence signals
|
---|
668 | cnt += arr->GetEntriesFast();
|
---|
669 |
|
---|
670 | // Remove all signals which are not in the valid digitization range
|
---|
671 | // (This is not the digitization window, but the region in which
|
---|
672 | // the analog channels contain usefull data)
|
---|
673 | // and which are shorter than the defined coincidence gate.
|
---|
674 | TIter Next(arr);
|
---|
675 | MDigitalSignal *ttl = 0;
|
---|
676 | while ((ttl=static_cast<MDigitalSignal*>(Next())))
|
---|
677 | {
|
---|
678 | if (ttl->GetLength()<fCoincidenceTime)
|
---|
679 | {
|
---|
680 | delete arr->Remove(ttl);
|
---|
681 | continue;
|
---|
682 | }
|
---|
683 |
|
---|
684 | if (ttl->GetStart()<min)
|
---|
685 | {
|
---|
686 | delete arr->Remove(ttl);
|
---|
687 | rmlo++;
|
---|
688 | continue;
|
---|
689 | }
|
---|
690 | if (ttl->GetStart()>max)
|
---|
691 | {
|
---|
692 | delete arr->Remove(ttl);
|
---|
693 | rmhi++;
|
---|
694 | continue;
|
---|
695 | }
|
---|
696 |
|
---|
697 | // Set trigger-channel index to keep this information
|
---|
698 | //ttl->SetIndex(j);
|
---|
699 | }
|
---|
700 |
|
---|
701 | // Remove the empty slots
|
---|
702 | arr->Compress();
|
---|
703 |
|
---|
704 | // cout << "ttls(j=" << j << "):";
|
---|
705 | // TObjArray *arr2 = static_cast<TObjArray*>(ttls[j]);
|
---|
706 | // TIter Nexty(arr);
|
---|
707 | // MDigitalSignal *ttly = 0;
|
---|
708 | // while ((ttly=static_cast<MDigitalSignal*>(Nexty())))
|
---|
709 | // {
|
---|
710 | // cout << "|"<< ttly->GetStart() << ", " << ttly->GetLength();
|
---|
711 | // }
|
---|
712 | // cout << endl;
|
---|
713 |
|
---|
714 |
|
---|
715 | // If we have at least one trigger keep the earliest one.
|
---|
716 | // FIXME: The triggers might be ordered in time automatically:
|
---|
717 | // To be checked!
|
---|
718 | // FIXME: Simulate trigger dead-time!
|
---|
719 | if (arr->GetEntriesFast()>0)
|
---|
720 | {
|
---|
721 | ttl = static_cast<MDigitalSignal*>(arr->RemoveAt(0));
|
---|
722 | // Set trigger-channel index to keep this information
|
---|
723 | ttl->SetIndex(j);
|
---|
724 | triggers.Add(ttl);
|
---|
725 | }
|
---|
726 |
|
---|
727 | // delete the allocated space
|
---|
728 | delete arr;
|
---|
729 | }
|
---|
730 |
|
---|
731 | // There are usually not enough entries that it is worth to search
|
---|
732 | // for the earliest instead of just sorting and taking the first one
|
---|
733 | // FIXME: This could be improved if checking for IsSortable
|
---|
734 | // is omitted
|
---|
735 | triggers.Sort();
|
---|
736 | // FIXME: Store triggers! (+ Reversed pixels?)
|
---|
737 |
|
---|
738 | SetTrigger(triggers.GetFirstTrigger(), triggers.GetFirstIndex());
|
---|
739 |
|
---|
740 | // No trigger issued. Go on.
|
---|
741 | if (triggers.GetNumSignals()==0)
|
---|
742 | {
|
---|
743 | if (rmlo>0 || rmhi>0)
|
---|
744 | *fLog << inf2 << GetNumExecutions() << ": " << rmlo << "/" << rmhi << " trigger out of valid range. No trigger raised." << endl;
|
---|
745 | return kTRUE;
|
---|
746 | }
|
---|
747 |
|
---|
748 | // Number of patches which have triggered out of the total number of
|
---|
749 | // Coincidence signals emitted. (If the total number is higher than
|
---|
750 | // the number of triggers either some triggers had to be removed or
|
---|
751 | // or a patch has emitted more than one trigger signal)
|
---|
752 | // FIXME: inf2?
|
---|
753 | *fLog << inf << GetNumExecutions() << ": ";
|
---|
754 | *fLog << setw(3) << triggers.GetNumSignals() << " triggers left out of ";
|
---|
755 | *fLog << setw(3) << cnt << " (" << rmlo << "/" << rmhi << " trigger out of valid range), T=" << fTrigger->GetVal();
|
---|
756 | *fLog << endl;
|
---|
757 |
|
---|
758 | //# Trigger characteristics: gate length (ns), min. overlapping time (ns),
|
---|
759 | //# amplitude and FWHM of (gaussian) single phe response for trigger:
|
---|
760 | //trigger_prop 3.0 0.25 1.0 2.0
|
---|
761 |
|
---|
762 | return kTRUE;
|
---|
763 | }
|
---|
764 |
|
---|
765 | Int_t MSimTrigger::PostProcess()
|
---|
766 | {
|
---|
767 | // patch_file.close();
|
---|
768 | // clipped_file.close();
|
---|
769 | // digital_file.close();
|
---|
770 | // ratescan_file.close();
|
---|
771 | return kTRUE;
|
---|
772 | }
|
---|
773 |
|
---|
774 |
|
---|
775 | // --------------------------------------------------------------------------
|
---|
776 | //
|
---|
777 | // FileNameRouteac: routeac.txt
|
---|
778 | // FileNameCoincidenceMap: coincidence.txt
|
---|
779 | // DiscriminatorTheshold: 3.5
|
---|
780 | // DigitalSignalLength: 8
|
---|
781 | // CoincidenceTime: 0.5
|
---|
782 | // SimulateElectronics: Yes
|
---|
783 | //
|
---|
784 | Int_t MSimTrigger::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
|
---|
785 | {
|
---|
786 | Bool_t rc = kFALSE;
|
---|
787 | if (IsEnvDefined(env, prefix, "FileNameRouteAC", print))
|
---|
788 | {
|
---|
789 | rc = kTRUE;
|
---|
790 | fNameRouteAC = GetEnvValue(env, prefix, "FileNameRouteAC", fNameRouteAC);
|
---|
791 | }
|
---|
792 |
|
---|
793 | if (IsEnvDefined(env, prefix, "FileNameCoincidenceMap", print))
|
---|
794 | {
|
---|
795 | rc = kTRUE;
|
---|
796 | fNameCoincidenceMap = GetEnvValue(env, prefix, "FileNameCoincidenceMap", fNameCoincidenceMap);
|
---|
797 | }
|
---|
798 |
|
---|
799 | if (IsEnvDefined(env, prefix, "DiscriminatorThreshold", print))
|
---|
800 | {
|
---|
801 | rc = kTRUE;
|
---|
802 | fDiscriminatorThreshold = GetEnvValue(env, prefix, "DiscriminatorThreshold", fDiscriminatorThreshold);
|
---|
803 | }
|
---|
804 |
|
---|
805 | if (IsEnvDefined(env, prefix, "DigitalSignalLength", print))
|
---|
806 | {
|
---|
807 | rc = kTRUE;
|
---|
808 | fDigitalSignalLength = GetEnvValue(env, prefix, "DigitalSignalLength", fDigitalSignalLength);
|
---|
809 | }
|
---|
810 |
|
---|
811 | if (IsEnvDefined(env, prefix, "CoincidenceTime", print))
|
---|
812 | {
|
---|
813 | rc = kTRUE;
|
---|
814 | fCoincidenceTime = GetEnvValue(env, prefix, "CoincidenceTime", fCoincidenceTime);
|
---|
815 | }
|
---|
816 |
|
---|
817 | if (IsEnvDefined(env, prefix, "SimulateElectronics", print))
|
---|
818 | {
|
---|
819 | rc = kTRUE;
|
---|
820 | fSimulateElectronics = GetEnvValue(env, prefix, "SimulateElectronics", fSimulateElectronics);
|
---|
821 | }
|
---|
822 |
|
---|
823 | if (IsEnvDefined(env, prefix, "MinMultiplicity", print))
|
---|
824 | {
|
---|
825 | rc = kTRUE;
|
---|
826 | fMinMultiplicity = GetEnvValue(env, prefix, "MinMultiplicity", fMinMultiplicity);
|
---|
827 | }
|
---|
828 |
|
---|
829 | if (IsEnvDefined(env, prefix, "CableDelay", print))
|
---|
830 | {
|
---|
831 | rc = kTRUE;
|
---|
832 | fCableDelay = GetEnvValue(env, prefix, "CableDelay", fCableDelay);
|
---|
833 | }
|
---|
834 |
|
---|
835 | if (IsEnvDefined(env, prefix, "CableDamping", print))
|
---|
836 | {
|
---|
837 | rc = kTRUE;
|
---|
838 | fCableDamping = GetEnvValue(env, prefix, "CableDamping", fCableDamping);
|
---|
839 | }
|
---|
840 |
|
---|
841 |
|
---|
842 |
|
---|
843 | return rc;
|
---|
844 | }
|
---|