| 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 | }
|
|---|