Changeset 9462 for trunk

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mbase/MLut.cc (modified) (6 diffs)
• mbase/MLut.h (modified) (2 diffs)
• mbase/MParList.cc (modified) (3 diffs)
• melectronics/MAvalanchePhotoDiode.cc (modified) (6 diffs)
• melectronics/MAvalanchePhotoDiode.h (modified) (3 diffs)
• mgeom/MGeomCamDwarf.cc (modified) (1 diff)
• mjobs/MJPedestal.cc (modified) (3 diffs)
• mjobs/MJSimulation.cc (modified) (1 diff)
• mpedestal/MPedCalcPedRun.cc (modified) (2 diffs)
• mraw/MRawRunHeader.cc (modified) (2 diffs)
• mraw/MRawRunHeader.h (modified) (1 diff)
• mreport/MReport.cc (modified) (2 diffs)
• mreport/MReportDrive.cc (modified) (1 diff)
• msimcamera/MSimAPD.cc (modified) (5 diffs)
• msimcamera/MSimTrigger.cc (modified) (19 diffs)
• msimcamera/MSimTrigger.h (modified) (1 diff)
• mtrigger/MTriggerPattern.h (modified) (1 diff)
• mtrigger/MTriggerPatternDecode.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2009/06/20 Thomas Bretz
+
+   * mbase/MLut.[h,cc]:
+     - renamed Default to DefaultCol
+     - added DefaultRow
+     - added Print
+     - added some const qualifiers
+
+   * mbase/MParList.cc:
+     - make sure a parlist is not added to itself
+     - improved some debug output
+
+   * melectronics/MAvalanchePhotoDiode.[h,cc]:
+     - keep a time-stamp fTime
+     - added functions to "relax/evolve" the apd over a given time
+
+   * mgeom/MGeomCamDwarf.cc:
+     - fixed a broken comment
+
+   * mjobs/MJPedestal.cc:
+     - improved output in case of error
+
+   * mjobs/MJSimulation.cc:
+     - implemented a PreCut (to gain processing time)
+     - implemented the correct type id in the file name
+
+   * mpedestal/MPedCalcPedRun.cc:
+     - check for the run-number only if real data
+
+   * mraw/MRawRunHeader.[h,cc]:
+     - added GetRuntypeChar
+     - fixed a bug which returned the runtype instead of the
+       telescopenumber
+
+   * mreport/MReport.cc:
+     - added a new report version number 200905170 (preliminray!)
+
+   * mreport/MReportDrive.cc:
+     - added a fix for the starguider sttaus in the report
+
+   * msimcamera/MSimAPD.cc:
+     - instead of a full initialization of the APD before each event
+       we now only simulate a time evolution which is needed to relax
+       a cell to one permille error. This is much faster.
+
+   * msimcamera/MSimTrigger.[h,cc]:
+     - allow switching off of the electronics trigger
+     - replaced the real coincidence by better and faster algorithm
+     - the new algorithm also supports multiplicity triggers
+     - therefore added fMinMultiplicity
+     - updated the handling and ouput in case of empty coincidence maps
+
+   * mtrigger/MTriggerPattern.h:
+     - added some information in the bit description
+
+   * mtrigger/MTriggerPatternDecode.cc:
+     - fixed the decoding for runs for which L1TPU was connected
+       instead of L1
+
+
 
  2009/06/18 Daniela Dorner

trunk/MagicSoft/Mars/mbase/MLut.cc

@@ -70 +70 @@
 //
 // Check if it is a default lut which would just map every entry to itself
-// An empty Lut is a default lut
-//
-Bool_t MLut::IsDefault() const
+// An empty Lut is a default column
+//
+Bool_t MLut::IsDefaultCol() const
 {
     if (IsEmpty())
@@ -92 +92 @@
 // Setup a default lut which just maps n-entris to themself
 //
-void MLut::SetDefault(UInt_t n)
+void MLut::SetDefaultCol(UInt_t n)
 {
     Delete();
@@ -106 +106 @@
     fMinEntries = 1;
     fMaxEntries = 1;
+
+    fMaxIndex = n;
+}
+
+// --------------------------------------------------------------------------
+//
+// Check if it is a default lut which would just map all entries to one.
+// An empty Lut is a default row
+//
+Bool_t MLut::IsDefaultRow() const
+{
+    if (IsEmpty())
+        return kTRUE;
+
+    if (GetEntriesFast()!=1)
+        return kFALSE;
+
+    const MArrayI &idx = GetRow(0);
+
+    // Loop over all rows
+    for (UInt_t x=0; x<idx.GetSize(); x++)
+        if (UInt_t(idx[x])!=x)
+            return kFALSE;
+
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Setup a default lut which maps all n-entris to one
+//
+void MLut::SetDefaultRow(UInt_t n)
+{
+    Delete();
+
+    MArrayI &idx = *new MArrayI(n);
+
+    for (UInt_t y=0; y<n; y++)
+        idx[y] = y;
+
+    Add(&idx);
+
+    fMinEntries = n;
+    fMaxEntries = n;
 
     fMaxIndex = n;
@@ -288 +333 @@
 // Write a lut to a stream.
 //
-Int_t MLut::WriteStream(ostream &out)
+Int_t MLut::WriteStream(ostream &out) const
 {
     const Int_t n = GetEntriesFast();
@@ -329 +374 @@
 // Write a lut to a file
 //
-Int_t MLut::WriteFile(const char *fname)
+Int_t MLut::WriteFile(const char *fname) const
 {
     TString expname(fname);
@@ -344 +389 @@
     return WriteStream(fout);
 }
+
+void MLut::Print(const Option_t *o) const
+{
+    gLog << all;
+    WriteStream(gLog);
+}

trunk/MagicSoft/Mars/mbase/MLut.h

@@ -36 +36 @@
     Bool_t HasConstantLength() const { return fMinEntries==fMaxEntries; }
     Bool_t IsEmpty() const { return fMaxEntries==0; }
-    Bool_t IsDefault() const;
+    Bool_t IsDefaultCol() const;
+    Bool_t IsDefaultRow() const;
 
     // MLut conversions
@@ -43 +44 @@
 
     // Setter
-    void SetDefault(UInt_t n);
+    void SetDefaultCol(UInt_t n);
+    void SetDefaultRow(UInt_t n);
 
     // MLut I/O
     Int_t ReadStream(istream &in);
-    Int_t WriteStream(ostream &out);
+    Int_t WriteStream(ostream &out) const;
 
     Int_t ReadFile(const char *fname);
-    Int_t WriteFile(const char *fname);
+    Int_t WriteFile(const char *fname) const;
+
+    // TObject
+    void Print(const Option_t *o="") const;
+    void Print(const Option_t *o, const Option_t *o2) const { Print(o); }
 
     ClassDef(MLut, 1) // A simple and fast easy-to-use look-up-table

trunk/MagicSoft/Mars/mbase/MParList.cc

@@ -184 +184 @@
         return kFALSE;
 
+    if (cont==this)
+    {
+        *fLog << err << dbginf << "Error: It is not allowed to add a parameter list to itself." << endl;
+        return kFALSE;
+    }
+
     //
     // Get Name of new container
@@ -223 +229 @@
         if (!fContainer->FindObject(where))
         {
-            *fLog << dbginf << "Error: Cannot find parameter container after which the new one should be added!" << endl;
+            *fLog << err << dbginf << "Error: Cannot find parameter container after which the new one should be added!" << endl;
             return kFALSE;
         }
@@ -230 +236 @@
     if (!cont->InheritsFrom(MParContainer::Class()))
     {
-        *fLog << dbginf << "Error: Cantainer MUST derive from MParContainer!" << endl;
+        *fLog << err << dbginf << "Error: Cantainer MUST derive from MParContainer!" << endl;
         return kFALSE;
     }

trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.cc

@@ -77 +77 @@
 APD::APD(Int_t n, Float_t prob, Float_t dt, Float_t rt)
     : fHist("APD", "", n, 0.5, n+0.5, n, 0.5, n+0.5),
-    fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt)
+    fCrosstalkProb(prob), fDeadTime(dt), fRecoveryTime(rt), fTime(-1)
 {
     fHist.SetDirectory(0);
+}
+
+// --------------------------------------------------------------------------
+//
+// This is the time the cell needs after a hit until less than threshold
+// (0.001 == one permille) of the signal is lost.
+//
+// If all cells of a G-APD have fired simultaneously and they are fired
+// once more after the relaxation time (with the defaultthreshold of 0.001)
+// the chip will show a signal which is one permille too small.
+//
+Float_t APD::GetRelaxationTime(Float_t threshold) const
+{
+    return fDeadTime-TMath::Log(threshold)*fRecoveryTime;
 }
 
@@ -197 +211 @@
 // The default time is 0.
 //
+// If you want t w.r.t. fTime use HitRandomCellRelative istead.
+//
 Float_t APD::HitRandomCell(Float_t t)
 {
@@ -219 +235 @@
 // If deadtime and recovery time are 0 then t-1 is set.
 //
+// Sets fTime to t
+//
 // The default time is 0.
 //
@@ -231 +249 @@
 
     fHist.SetEntries(1);
+
+    fTime = t;
 }
 
@@ -238 +258 @@
 // by calling HitCellImp with time t-gRandom->Exp(n/rate) with n being
 // the total number of cells.
+//
+// Sets fTime to t
+//
 // The default time is 0.
 //
@@ -253 +276 @@
     for (int x=1; x<=nx; x++)
         for (int y=1; y<=ny; y++)
-        {
             HitCellImp(x, y, t-MMath::RndmExp(f));
-        }
+
+    fTime = t;
+}
+
+// --------------------------------------------------------------------------
+//
+// Evolve the chip from fTime to fTime+dt if it with a given frequency
+// freq. Returns the total signal "recorded".
+//
+// fTime is set to the fTime+dt
+//
+// If you want to evolve over a default relaxation time (relax the chip
+// from a signal) use Relax instead.
+//
+Float_t APD::Evolve(Double_t freq, Double_t dt)
+{
+    const Double_t avglen = 1./freq;
+
+    const Double_t end  = fTime+dt;
+
+    Float_t hits = 0;
+
+    Double_t time = fTime;
+    while (1)
+    {
+        time += MMath::RndmExp(avglen);
+        if (time>end)
+            break;
+
+        hits += HitRandomCell(time);
+    }
+
+    fTime = end;
+
+    return hits;
 }
 

trunk/MagicSoft/Mars/melectronics/MAvalanchePhotoDiode.h

@@ -15 +15 @@
     Float_t fRecoveryTime;   // Recoverytime after Deadtime (1-exp(-t/fRecoveryTime)
 
+    Float_t fTime;           // A user settable time of the system
+
     Float_t HitCellImp(Int_t x, Int_t y, Float_t t=0);
 
@@ -22 +24 @@
     Float_t HitCell(Int_t x, Int_t y, Float_t t=0);
     Float_t HitRandomCell(Float_t t=0);
+    Float_t HitRandomCellRelative(Float_t t=0) { return HitRandomCell(fTime+t); }
 
     void FillEmpty(Float_t t=0);
     void FillRandom(Float_t rate, Float_t t=0);
+
+    void Init(Float_t rate) { if (fTime<0) FillRandom(rate); else Relax(rate); }
 
     Int_t CountDeadCells(Float_t t=0) const;
@@ -35 +40 @@
     Float_t GetDeadTime() const { return fDeadTime; }
     Float_t GetRecoveryTime() const { return fRecoveryTime; }
+    Float_t GetTime() const { return fTime; }
+
+    Float_t GetRelaxationTime(Float_t threshold=0.001) const;
+
+    Float_t GetLastHit() const { return fHist.GetMaximum(); }
+
+    void SetTime(Float_t tm) { fTime=tm; }
+    void IncreaseTime(Float_t dt) { fTime += dt; }
+
+    Float_t Evolve(Double_t freq, Double_t dt);
+    Float_t Relax(Double_t freq, Float_t threshold=0.001) { return Evolve(freq, GetRelaxationTime(threshold)); }
 
     void Draw(Option_t *o="") { fHist.Draw(o); }

trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.cc

@@ -38 +38 @@
 // hexagonal camera, but the edges filled with additional pixels
 // inside a circle.
-//    MGeomCamDwarf cam(
+//    MGeomCamDwarf cam(209.5, 13.2);
 //
 ////////////////////////////////////////////////////////////////////////////

trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

@@ -1154 +1154 @@
     if (taskenv.GetNumExecutions()<fMinEvents)
     {
-        *fLog << err << GetDescriptor() << ": Failed. Less than the required " << fMinEvents << " evts processed." << endl;
+        *fLog << err << GetDescriptor() << ": Failed. Less (" << taskenv.GetNumExecutions() << ") than the required " << fMinEvents << " pedestal events extracted." << endl;
         return kFALSE;
     }
@@ -1160 +1160 @@
     if (fIsPulsePosCheck && pulcam.GetNumEvents()<fMinCosmics)
     {
-        *fLog << err << GetDescriptor() << ": Failed. Less than the required " << fMinCosmics << " cosmics evts processed." << endl;
+        *fLog << err << GetDescriptor() << ": Failed. Less (" << pulcam.GetNumEvents() << ") than the required " << fMinCosmics << " cosmics evts processed." << endl;
         return kFALSE;
     }
@@ -1166 +1166 @@
     if (fPedestalCamOut.GetNumEvents()<fMinPedestals)
     {
-        *fLog << err << GetDescriptor() << ": Failed. Less than the required " << fMinPedestals << " pedetsl evts processed." << endl;
+        *fLog << err << GetDescriptor() << ": Failed. Less (" << fPedestalCamOut.GetNumEvents() << ") than the required " << fMinPedestals << " pedestal evts processed." << endl;
         return kFALSE;
     }

trunk/MagicSoft/Mars/mjobs/MJSimulation.cc

@@ -551 +551 @@
     MSimCalibrationSignal simcal;
     simcal.SetNameGeomCam("GeomCones");
+
+    // Sky Quality Meter: 21.82M = 2.02e-4cd/m^2
+    // 1cd = 12.566 lm / 4pi sr
 
     // FIXME: Simulate photons before cones and QE!

trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc

@@ -1 +1 @@
 /* ======================================================================== *\
-! $Name: not supported by cvs2svn $:$Id: MPedCalcPedRun.cc,v 1.54 2008-11-12 16:04:18 tbretz Exp $
+! $Name: not supported by cvs2svn $:$Id: MPedCalcPedRun.cc,v 1.55 2009-06-20 09:14:33 tbretz Exp $
 ! --------------------------------------------------------------------------
 !
@@ -320 +320 @@
 //
 Bool_t MPedCalcPedRun::IsPedBitSet()
-{
-    if (fRunHeader->GetTelescopeNumber()==1 && fRunHeader->GetRunNumber()<gkFirstRunWithFinalBits)
+{     
+    if (!fRunHeader->IsMonteCarloRun() && fRunHeader->GetTelescopeNumber()==1 && fRunHeader->GetRunNumber()<gkFirstRunWithFinalBits)
         return kFALSE;
 

trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc

@@ -940 +940 @@
 // to print it as readable text.
 //
-const char *MRawRunHeader::GetRunTypeStr() const
+const Char_t *MRawRunHeader::GetRunTypeStr() const
 {
     switch (fRunType)
@@ -960 +960 @@
     default:
         return "<unknown>";
+    }
+}
+
+const Char_t MRawRunHeader::GetRunTypeChar() const
+{
+    switch (fRunType&0xff)
+    {
+    case kRTData:
+    case kRTPointRun:
+        return 'D';
+    case kRTPedestal:
+        return 'P';
+    case kRTCalibration:
+        return 'C';
+    case kRTLinearity:
+        return 'N';
+    default:
+        return ' ';
     }
 }

trunk/MagicSoft/Mars/mraw/MRawRunHeader.h

@@ -126 +126 @@
     UInt_t   GetFileID() const            { return fRunNumber>1000000?(fRunNumber%1000000)*1000+(fFileNumber%1000):fRunNumber; }
     TString  GetStringID() const;
-    UShort_t GetTelescopeNumber() const   { return fRunType; }
+    UShort_t GetTelescopeNumber() const   { return fTelescopeNumber; }
     UShort_t GetRunType() const           { return fRunType; }
     const Char_t *GetRunTypeStr() const;
+    const Char_t  GetRunTypeChar() const;
     const Char_t *GetProjectName() const  { return fProjectName; }
     const Char_t *GetSourceName() const   { return fSourceName; }

trunk/MagicSoft/Mars/mreport/MReport.cc

@@ -169 +169 @@
 //    200805190  | 54711.5 |         |  200809030
 //    200812040  | 54814.5 |         |  200812140
+//    200902210  | 54968.5 |         |  200905170
 //
 Int_t MReport::Interprete(TString &str, const MTime &start, const MTime &stop, Int_t ver)
@@ -208 +209 @@
         ver=200812140;
 
+    if (ver==200902210 && GetMjd()>54968.5)
+        ver=200905170;
+
     // Interprete body (contents) of report
     const Int_t rc = InterpreteBody(str, ver);

trunk/MagicSoft/Mars/mreport/MReportDrive.cc

@@ -135 +135 @@
     }
 
+    if (ver>=200805170)
+    {
+        Int_t dummy; // Starguider switched on or not
+        n=sscanf(str.Data(), "%d %n", &dummy, &len);
+        if (n!=1)
+        {
+            *fLog << warn << "WARNING - Not enough arguments." << endl;
+            return kCONTINUE;
+        }
+
+        str.Remove(0, len);
+        str = str.Strip(TString::kBoth);
+    }
+
     return str.IsNull() ? kTRUE : kCONTINUE;
 }

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

trunk/MagicSoft/Mars/mtrigger/MTriggerPattern.h

@@ -12 +12 @@
 public:
     enum Pattern_t {
-        kTriggerLvl1 = BIT(0), //  1:
-        kCalibration = BIT(1), //  2: Pulse Trigger
-        kTriggerLvl2 = BIT(2), //  4: LUT Pseudo Size selection
-        kPedestal    = BIT(3), //  8:
-        kPinDiode    = BIT(4), // 16:
-        kSumTrigger  = BIT(5)  // 32: Flag for an event taken with sum trigger
-        //kUnused      = BIT(6)
-        //kUnused      = BIT(7)
+//        kUndefined1  = BIT(0), //  1   1: Level 1 from L2 board
+        kTriggerLvl1 = BIT(0), //  1   1: Level 1 from L2 board
+        kCalibration = BIT(1), //  2   2: Pulse Trigger
+        kTriggerLvl2 = BIT(2), //  4   4: LUT Pseudo Size selection
+        kPedestal    = BIT(3), //  8   8: Artificial pedestal event
+        kPinDiode    = BIT(4), // 10  16:
+        kSumTrigger  = BIT(5), // 20  32: Flag for an event taken with sum trigger
+//        kTriggerLvl1 = BIT(6), // 40  64: Trigger lvl1 directly from L1 without going through L2
+        kUndefined1  = BIT(6), // 40  64: Trigger lvl1 directly from L1 without going through L2
+        kUndefined2  = BIT(7)  // 80 128: Undefined?
     };
 

trunk/MagicSoft/Mars/mtrigger/MTriggerPatternDecode.cc

@@ -106 +106 @@
         return kTRUE;
 
-    UInt_t pattern = ~fEvtHeader->GetTriggerID();
+    const UInt_t pattern = ~fEvtHeader->GetTriggerID();
 
     // The the trigger pattern is currently written with inverted bits,
@@ -133 +133 @@
     fPattern->fUnprescaled = (pattern>>8) & 0xff;
 
+    // This is a workaround for the new scheme in which L1TPU (the signal
+    // comming directly from the L1 is connected, but the L1 (routed
+    // over L2 is disconnected)
+    if (!fRunHeader->IsMonteCarloRun() && fRunHeader->GetTelescopeNumber()==1 &&
+        fRunHeader->GetRunNumber()>1006246)
+    {
+        fPattern->fPrescaled   |= (pattern>> 6)&1;
+        fPattern->fUnprescaled |= (pattern>>14)&1;
+    }
+
     return kTRUE;
 }