Changeset 9462 for trunk/MagicSoft/Mars/msimcamera

Timestamp:

06/20/09 10:14:33 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/msimcamera

Files:

• MSimAPD.cc (modified) (5 diffs)
• MSimTrigger.cc (modified) (19 diffs)
• MSimTrigger.h (modified) (1 diff)

trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc

@@ -206 +206 @@
         return kERROR;
     }
-
+/*
     for (UInt_t idx=0; idx<npix; idx++)
     {
@@ -212 +212 @@
         static_cast<APD*>(fAPDs.UncheckedAt(idx))->FillRandom(freq, fStat->GetTimeFirst());
     }
+*/
+
+    // This tries to initialize dead and relaxing cells properly. If
+    // the APD has not been initialized before the chip is randomsly
+    // filled, otherwise a time window of the default relaxing time
+    // is simulated, so that the previous influence is less than a permille.
+    for (UInt_t idx=0; idx<npix; idx++)
+    {
+        const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
+        static_cast<APD*>(fAPDs.UncheckedAt(idx))->Init(freq);
+    }
 
     // Get number of photons
@@ -223 +234 @@
 
         // Get arrival time of photon and idx
-        const Double_t t = ph.GetTime();
+        const Double_t t = ph.GetTime()-fStat->GetTimeFirst();
         const Int_t  idx = ph.GetTag();
         if (idx<0)
@@ -239 +250 @@
         // Simulate hitting the APD (the signal height in effective
         // "number of photons" is returned)
-        const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCell(t);
+        const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCellRelative(t);
 
         // FIXME: Make a proper simulation of the excess noise!!!
@@ -248 +259 @@
     }
 
+    // Now we have to shift the evolved time of all APDs to the end of our
+    // simulated time.
+    for (UInt_t idx=0; idx<npix; idx++)
+        static_cast<APD*>(fAPDs.UncheckedAt(idx))->IncreaseTime(fStat->GetTimeLast());
+
     return kTRUE;
 }

trunk/MagicSoft/Mars/msimcamera/MSimTrigger.cc

@@ -42 +42 @@
 //
 // With a second look-up table fCoincidenceMap the analog channels are
-// checked for coincidences. The coincidense must at least be of the length
+// checked for coincidences. The coincidence must at least be of the length
 // defined by fCoincidenceTime. The earliest coincide is then stored as
 // trigger position.
+//
+// If a minimum multiplicity m is given, m signals above threshold
+// in the coincidence patterns are enough to emit a trigger signal.
 //
 //
@@ -97 +100 @@
     fEvtHeader(0), fElectronicNoise(0), fGain(0),
     fDiscriminatorThreshold(-1), fDigitalSignalLength(8), fCoincidenceTime(0.5),
-    fShiftBaseline(kTRUE), fUngainSignal(kTRUE)
+    fShiftBaseline(kTRUE), fUngainSignal(kTRUE), fSimulateElectronics(kTRUE),
+    fMinMultiplicity(-1)
 {
     fName  = name  ? name  : "MSimTrigger";
@@ -107 +111 @@
 // Take two TObjArrays with a collection of digital signals.
 // Every signal from one array is compared with any from the other array.
-// For all signals whihc overlaps and which have an overlap time >gate
+// For all signals which overlap and which have an overlap time >gate
 // a new digital signal is created storing start time and length of overlap.
 // They are collected in a newly allocated TObjArray. A pointer to this array
 // is returned.
 //
-// Th euser gains owenership of the object, ie.e., the user is responsible of
+// The user gains owenership of the object, i.e., the user is responsible of
 // deleting the memory.
 //
-TObjArray *MSimTrigger::CalcCoincidence(const TObjArray &arr1, const TObjArray &arr2, Float_t gate) const
+TObjArray *MSimTrigger::CalcCoincidence(const TObjArray &arr1, const TObjArray &arr2/*, Float_t gate*/) const
 {
     TObjArray *res = new TObjArray;
@@ -131 +135 @@
         {
             MDigitalSignal *ttl = new MDigitalSignal(*ttl1, *ttl2);
-
+            /*
             if (ttl->GetLength()<=gate)
             {
@@ -137 +141 @@
                 continue;
             }
-
+            */
             res->Add(ttl);
         }
@@ -143 +147 @@
 
     res->SetOwner();
+
+    return res;
+}
+
+class Edge : public TObject
+{
+private:
+    Double_t fEdge;
+    Int_t    fRising;
+
+public:
+    Edge(Double_t t, Int_t rising) : fEdge(t), fRising(rising) { }
+    Bool_t IsSortable() const { return kTRUE; }
+    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; }
+
+    Int_t IsRising() const { return fRising; }
+    Double_t GetEdge() const { return fEdge; }
+};
+
+// --------------------------------------------------------------------------
+//
+// Calculate a multiplicity trigger on the given array(s). The idx-array
+// conatins all channels which should be checked for coincidences
+// and the ttls array conatins the arrays with the digital signals.
+//
+// For the windows in which more or euqal than threshold channels have
+// a high signal a new MDigitalSignal is created.  newly allocated
+// array with a collection of these trigger signals is returned.
+//
+TObjArray *MSimTrigger::CalcMinMultiplicity(const MArrayI &idx, const TObjArray &ttls, Int_t threshold) const
+{
+    // Create a new array for the rising and falling edges of the signals
+    TObjArray times;
+    times.SetOwner();
+
+    // Fill the array with edges from all digital signals of all our channels
+    for (UInt_t k=0; k<idx.GetSize(); k++)
+    {
+        TObjArray *arr = static_cast<TObjArray*>(ttls[idx[k]]);
+
+        TIter Next(arr);
+        MDigitalSignal *ttl = 0;
+        while ((ttl=static_cast<MDigitalSignal*>(Next())))
+        {
+            times.Add(new Edge(ttl->GetStart(),  1));
+            times.Add(new Edge(ttl->GetEnd(),   -1));
+        }
+    }
+
+    // Sort them in time
+    times.Sort();
+
+    // Start with no channel active
+    Int_t lvl = 0;
+
+    TObjArray *res = new TObjArray;
+    res->SetOwner();
+
+    // First remove all edges which do not change the status
+    // "below threshold" or "above threshold"
+    for (int i=0; i<times.GetEntriesFast(); i++)
+    {
+        // Get i-th edge
+        const Edge &e = *static_cast<Edge*>(times.UncheckedAt(i));
+
+        // Claculate what the number of active channels after the edge is
+        const Int_t lvl1 = lvl + e.IsRising();
+
+        // Remove edge if number of active channels before and after the
+        // edge lower is lower than the threshold or higher than
+        // the threshold
+        if (lvl+1<threshold || lvl-1>=threshold)
+            delete times.RemoveAt(i);
+
+        // keep the (now) "previous" level
+        lvl = lvl1<0 ? 0 : lvl1;
+    }
+
+    // Remove the empty slots from the array
+    times.Compress();
+
+    //
+    for (int i=0; i<times.GetEntriesFast()-1; i++)
+    {
+        // get the current edge
+        const Edge &e0 = *static_cast<Edge*>(times.UncheckedAt(i));
+
+        // go ahead if this is a falling edge
+        if (e0.IsRising()!=1)
+            continue;
+
+        // get the following edge (must be a falling edge now)
+        const Edge &e1 = *static_cast<Edge*>(times.UncheckedAt(i+1));
+
+        // calculate the length of the digital signal
+        const Double_t len = e1.GetEdge()-e0.GetEdge();
+
+        // Create a digital trigger signal
+        MDigitalSignal *ds = new MDigitalSignal(e0.GetEdge(), len);
+        //ds->SetIndex(lvl);
+        res->Add(ds);
+    }
 
     return res;
@@ -231 +337 @@
     }
 
-    if (fElectronicNoise && !fRouteAC.IsEmpty() && !fRouteAC.IsDefault())
+    if (fElectronicNoise && !fRouteAC.IsEmpty() && !fRouteAC.IsDefaultCol())
     {
         // FIXME: Apply to analog channels when summing
@@ -238 +344 @@
     }
 
-    if (fGain && !fRouteAC.IsEmpty() && !fRouteAC.IsDefault())
+    if (fGain && !fRouteAC.IsEmpty() && !fRouteAC.IsDefaultCol())
     {
         // FIXME: Apply to analog channels when summing
@@ -255 +361 @@
         *fLog << "Using " << fNameRouteAC << " for re-routing/summing of analog channels before discriminator." << endl;
 
-    if (fCoincidenceMap.IsEmpty())
+    if (fCoincidenceMap.IsEmpty() && fMinMultiplicity==0)
         *fLog << "No coincidences of digital channels will be checked. Signal-above-threshold trigger applied." << endl;
     else
-        *fLog << "Using " << fNameCoincidenceMap << " to check for coincidences of the digital channels." << endl;
+    {
+        *fLog << "Using ";
+        if (fCoincidenceMap.IsEmpty())
+            *fLog << "the whole camera";
+        else
+            *fLog << "patterns from " << fNameCoincidenceMap;
+        *fLog << " to check for ";
+        if (fMinMultiplicity==0)
+            *fLog << "coincidences of the digital channels." << endl;
+        else
+            *fLog << fMinMultiplicity << " multiplicity." << endl;
+        }
 
     *fLog << "Using discriminator threshold of " << fDiscriminatorThreshold << endl;
@@ -265 +382 @@
 }
 
+/*
+class MDigitalChannel : public TObjArray
+{
+private:
+    TObjArray fArray;
+
+public:
+    MDigitalSignal *GetSignal(UInt_t i) { return static_cast<MDigitalSignal*>(fArray[i]); }
+
+};
+*/
+
+#include "MCamEvent.h"
+class MTriggerSignal : public MParContainer, public MCamEvent
+{
+private:
+    TObjArray fSignals;
+
+public:
+    MTriggerSignal() { fSignals.SetOwner(); }
+
+    void Add(MDigitalSignal *signal) { fSignals.Add(signal); }
+
+    MDigitalSignal *GetSignal(UInt_t i) { return static_cast<MDigitalSignal*>(fSignals[i]); }
+
+    void Sort() { fSignals.Sort(); }
+
+    Int_t GetNumSignals() const { return fSignals.GetEntriesFast(); }
+
+    Float_t GetFirstTrigger() const
+    {
+        MDigitalSignal *sig = static_cast<MDigitalSignal*>(fSignals[0]);
+        return sig ? sig->GetStart() : -1;
+    }
+    Bool_t GetPixelContent(Double_t&, Int_t, const MGeomCam&, Int_t) const
+    {
+        switch (1)
+        {
+        case 1: // yes/no
+        case 2: // first time
+        case 3: // length
+        case 4: // n
+            break;
+        }
+
+        return kTRUE;
+    }
+    void DrawPixelContent(Int_t) const
+    {
+    }
+};
+
+
+void MSimTrigger::SetTrigger(Double_t pos)
+{
+    // FIXME: Jitter! (Own class?)
+    fTrigger->SetVal(pos);
+    fTrigger->SetReadyToSave();
+
+    // Trigger pattern to be set
+    // FIXME: Make flexible
+    const UInt_t pat = (UInt_t)~(MTriggerPattern::kTriggerLvl1 | (MTriggerPattern::kTriggerLvl1<<8));
+
+    // Flag this event as triggered by the lvl1 trigger (FIXME?)
+    fEvtHeader->SetTriggerPattern(pos<0 ? 0 : pat);
+    fEvtHeader->SetReadyToSave();
+}
+
 // --------------------------------------------------------------------------
 //
@@ -270 +455 @@
 {
     // Invalidate trigger
-    fTrigger->SetVal(-1);
+    //fTrigger->SetVal(-1);
+    // Calculate the minimum and maximum time for a valid trigger
+    const Double_t freq    = fRunHeader->GetFreqSampling()/1000.;
+    const Float_t  nsamp   = fRunHeader->GetNumSamplesHiGain();
+    const Float_t  pulspos = fPulsePos->GetVal()/freq;
+
+    // Valid range in units of bins
+    const Float_t min = fCamera->GetValidRangeMin()+pulspos;
+    const Float_t max = fCamera->GetValidRangeMax()-(nsamp-pulspos);
+
+    if (!fSimulateElectronics)
+    {
+        SetTrigger(min);
+        return kTRUE;
+    }
 
     // ================== Simulate channel bundling ====================
@@ -297 +496 @@
             {
                 const UInt_t idx = row[j];
-                (*patches)[i].AddSignal((*fCamera)[idx]);
+
+                // FIXME: Shrinking the mapping table earlier (e.g.
+                //        ReInit) would avoid a lot of if's
+                if (idx<fCamera->GetNumChannels())
+                    (*patches)[i].AddSignal((*fCamera)[idx]);
             }
         }
@@ -329 +532 @@
     //        omitted in the loop
     if (fCoincidenceMap.IsEmpty())
-        fCoincidenceMap.SetDefault(npatch);
-
-    // Calculate the minimum and maximum time for a valid trigger
-    const Double_t freq    = fRunHeader->GetFreqSampling()/1000.;
-    const Float_t  nsamp   = fRunHeader->GetNumSamplesHiGain();
-    const Float_t  pulspos = fPulsePos->GetVal()/freq;
-
-    // Valid range in units of bins
-    const Float_t min = fCamera->GetValidRangeMin()+pulspos;
-    const Float_t max = fCamera->GetValidRangeMax()-(nsamp-pulspos);
+    {
+        if (fMinMultiplicity>0)
+            fCoincidenceMap.SetDefaultRow(npatch);
+        else
+            fCoincidenceMap.SetDefaultCol(npatch);
+    }
 
     // Create an array for the individual triggers
-    TObjArray triggers;
-    triggers.SetOwner();
+    MTriggerSignal triggers;
 
     Int_t cnt  = 0;
@@ -352 +550 @@
         const MArrayI &idx = fCoincidenceMap.GetRow(j);
 
-        // Start with a copy of the first coincidence channel
-        TObjArray *arr = static_cast<TObjArray*>(ttls[idx[0]]->Clone());
-        arr->SetOwner();
-
-        // compare to all other channels in this coincidence patch, one by one
-        for (UInt_t k=1; k<idx.GetSize() && arr->GetEntriesFast()>0; k++)
-        {
-            TObjArray *res = CalcCoincidence(*arr, *static_cast<TObjArray*>(ttls[idx[k]]),
-                                             fCoincidenceTime);
-
-            // Delete the original array and keep the new one
-            delete arr;
-            arr = res;
-        }
+        TObjArray *arr = 0;
+
+        if (fMinMultiplicity>0)
+        {
+            arr = CalcMinMultiplicity(idx, ttls, fMinMultiplicity);
+        }
+        else
+        {
+            arr = CalcMinMultiplicity(idx, ttls, idx.GetSize());
+            /*
+            // Start with a copy of the first coincidence channel
+            arr = static_cast<TObjArray*>(ttls[idx[0]]->Clone());
+            arr->SetOwner();
+
+            // compare to all other channels in this coincidence patch, one by one
+            for (UInt_t k=1; k<idx.GetSize() && arr->GetEntriesFast()>0; k++)
+            {
+                TObjArray *res = CalcCoincidence(*arr, *static_cast<TObjArray*>(ttls[idx[k]]));//, fCoincidenceTime);
+
+                // Delete the original array and keep the new one
+                delete arr;
+                arr = res;
+            }*/
+        }
+
+        // Count the number of totally emitted coincidence signals
+        cnt += arr->GetEntriesFast();
 
         // Remove all signals which are not in the valid digitization range
         // (This is not the digitization window, but the region in which
         //  the analog channels contain usefull data)
+        // and which are shorter than the defined coincidence gate.
         TIter Next(arr);
         MDigitalSignal *ttl = 0;
         while ((ttl=static_cast<MDigitalSignal*>(Next())))
         {
+            if (ttl->GetLength()<fCoincidenceTime)
+            {
+                delete arr->Remove(ttl);
+                continue;
+            }
+
             if (ttl->GetStart()<min)
             {
                 delete arr->Remove(ttl);
                 rmlo++;
+                continue;
             }
             if (ttl->GetStart()>max)
@@ -383 +602 @@
                 delete arr->Remove(ttl);
                 rmhi++;
+                continue;
             }
+
+            // Set trigger channel index
+            ttl->SetIndex(j);
         }
 
@@ -389 +612 @@
         arr->Compress();
 
-        cnt += arr->GetEntriesFast();
-
         // If we have at least one trigger keep the earliest one.
-        // FIXME: The triggers should be ordered in time automatically: To be checked!
+        // FIXME: The triggers might be ordered in time automatically:
+        //        To be checked!
         // FIXME: Simulate trigger dead-time!
         if (arr->GetEntriesFast()>0)
-            triggers.Add(arr->RemoveAt(0));
+            triggers.Add(static_cast<MDigitalSignal*>(arr->RemoveAt(0)));
 
         // delete the allocated space
         delete arr;
-    }
-
-    // No trigger issued. Go on.
-    if (triggers.GetEntriesFast()==0)
-    {
-        if (rmlo>0 || rmhi>0)
-            *fLog << inf2 << rmlo << "/" << rmhi << " trigger out of valid range. No trigger raised." << endl;
-        return kTRUE;
     }
 
@@ -414 +628 @@
     //         is omitted
     triggers.Sort();
-
-    // FIXME: Jitter! (Own class?)
-    fTrigger->SetVal(static_cast<MDigitalSignal*>(triggers[0])->GetStart());
-    fTrigger->SetReadyToSave();
-
-    // Flag this event as triggered by the lvl1 trigger (FIXME?)
-    fEvtHeader->SetTriggerPattern((UInt_t)~(MTriggerPattern::kTriggerLvl1 | (MTriggerPattern::kTriggerLvl1<<8)));
-    fEvtHeader->SetReadyToSave();
-
-    // inf2?
-    *fLog << inf;
-    *fLog << cnt << " triggers left in " << triggers.GetEntriesFast() << " patches (" << rmlo << "/" << rmhi << " trigger out of valid range), T=" << fTrigger->GetVal();
+    // FIXME: Store triggers! (+ Reversed pixels?)
+
+    SetTrigger(triggers.GetFirstTrigger());
+
+    // No trigger issued. Go on.
+    if (triggers.GetNumSignals()==0)
+    {
+        if (rmlo>0 || rmhi>0)
+            *fLog << inf2 << GetNumExecutions() << ": " << rmlo << "/" << rmhi << " trigger out of valid range. No trigger raised." << endl;
+        return kTRUE;
+    }
+
+    // Number of patches which have triggered out of the total number of
+    // Coincidence signals emitted. (If the total number is higher than
+    // the number of triggers either some triggers had to be removed or
+    // or a patch has emitted more than one trigger signal)
+    // FIXME: inf2?
+    *fLog << inf << GetNumExecutions() << ": ";
+    *fLog << setw(3) << triggers.GetNumSignals() << " triggers left out of ";
+    *fLog << setw(3) << cnt << " (" << rmlo << "/" << rmhi << " trigger out of valid range), T=" << fTrigger->GetVal();
     *fLog << endl;
+
+    //# Trigger characteristics: gate length (ns), min. overlapping time (ns),
+    //# amplitude and FWHM of (gaussian) single phe response for trigger:
+    //trigger_prop 3.0 0.25 1.0 2.0
 
     return kTRUE;
@@ -438 +664 @@
 // DigitalSignalLength:     8
 // CoincidenceTime:         0.5
+// SimulateElectronics:     Yes
 //
 Int_t MSimTrigger::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
@@ -472 +699 @@
     }
 
+    if (IsEnvDefined(env, prefix, "SimulateElectronics", print))
+    {
+        rc = kTRUE;
+        fSimulateElectronics = GetEnvValue(env, prefix, "SimulateElectronics", fSimulateElectronics);
+    }
+
+    if (IsEnvDefined(env, prefix, "MinMultiplicity", print))
+    {
+        rc = kTRUE;
+        fMinMultiplicity = GetEnvValue(env, prefix, "MinMultiplicity", fMinMultiplicity);
+    }
+
     return rc;
 }

trunk/MagicSoft/Mars/msimcamera/MSimTrigger.h

@@ -40 +40 @@
     Bool_t  fShiftBaseline;             // Shift the baseline back to 0 for the threshold (needs ElectronicNoise [MPedestalCam])
     Bool_t  fUngainSignal;              // "Remove" the gain from the signal (needs Gain [MPedestalCam])
+    Bool_t  fSimulateElectronics;       // If the electronics is not simulated the trigger is set artificially to the first photon arrived
+
+    Int_t   fMinMultiplicity;           // N out of M
 
     // MSimTrigger
-    TObjArray *CalcCoincidence(const TObjArray &arr1, const TObjArray &arr2, Float_t gate=0) const;
+    TObjArray *CalcCoincidence(const TObjArray &arr1, const TObjArray &arr2/*, Float_t gate=0*/) const;
+    TObjArray *CalcMinMultiplicity(const MArrayI &idx, const TObjArray &ttls, Int_t threshold) const;
+    TObjArray *CalcCoincidences(const MArrayI &idx, const TObjArray &ttls) const;
+    void SetTrigger(Double_t pos);
 
     // MTask