source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/electronics/src/ChipTrackingTrgEngine0.cc @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 20.7 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: ChipTrackingTrgEngine0.cc 2920 2011-06-11 06:46:41Z mabl $
3// R. Pesce created
4
5
6//_____________________________________________________________________________
7//
8//  Chip Tracking Trigger Engine 0
9//  ==============================
10//
11//  Tracking trigger built in front end chip.
12//  The algorithm searches for sequences of pixels (over the threshold)
13//  consecutive and adjacent. If the track length is over a threshold,
14//  trigger occurs.
15//
16//  In the following figure there is explained the trigger concept:
17//
18//Begin_Html
19/*
20  <img src="./../pics/trigger_concept.gif">
21*/
22//End_Html
23//
24//  The algorithm searches also between track segments of two nearby chips:
25//  if two track segments are consecutive in time and the total length is
26//  over a threshold (different from the previous one), trigger occurs
27//
28//  Parameters:
29//
30//  fThreshold: threshold on pixel hits in a gtu
31//  fMinTrackLength: min length of the valid tracks
32//  fMaxTrackLength: max length of the valid tracks
33//  fMinTriggerTrackLength: min track length to have a trigger (single chip)
34//  fMinTriggerTwoLength: min total track length in 2 nearby chips
35//
36//  fAcceptHole [bool] : an hole in a single track segment is or not accepted
37//  fOnlyWithSignal [bool] : processes or not only chips with at least a signal
38//
39//  Special class for studying background fake rate trigger:
40//  The algorithm counts the number of tracks of the various trigger lengths,
41//  fill a specific trigger word and dump the numbers of the recognized tracks
42//
43
44#include "ChipTrackingTrgEngine0.hh"
45#include "MacroCellData.hh"
46#include "MacroCell.hh"
47#include "ElementaryCell.hh"
48#include "FrontEndChip.hh"
49#include "BoolAlgebra.hh"
50#include "ChipTrackSegment.hh"
51#include "EusoElectronics.hh"
52#include "Photomultiplier.hh"
53#include "EusoDetector.hh"
54#include "TMath.h"
55#include "ERunParameters.hh"
56#include "EChipTrackTriggerDataAdder.hh"
57#include "EChipTriggParsFiller.hh"
58#include "Config.hh"
59
60using namespace TMath;
61
62// trigger pattern types
63const Int_t kTrigPattern[32] = {BIT(0),BIT(1),BIT(2),BIT(3),BIT(4),BIT(5),BIT(6),BIT(7),
64                                BIT(8),BIT(9),BIT(10),BIT(11),BIT(12),BIT(13),BIT(14),BIT(15),
65                                BIT(16),BIT(17),BIT(18),BIT(19),BIT(20),BIT(21),BIT(22),BIT(23),
66                                BIT(24),BIT(25),BIT(26),BIT(27),BIT(28),BIT(29),BIT(30),BIT(31)};
67
68ClassImp(ChipTrackingTrgEngine0)
69
70//______________________________________________________________________________
71ChipTrackingTrgEngine0::ChipTrackingTrgEngine0() :
72TriggerEngine(string("ChipTrackingTrgEngine0"), kChipTrackingTrigger0 ) {
73    //
74    // Constructor
75    //
76
77    fMinTrackLength = 2;
78    fMaxTrackLength = 16;
79    fMinTriggerTrackLength = 2;
80    fMinTriggerTwoLength = 5;
81    fMaxTwoLength = 18;
82
83    fAcceptHole = kTRUE;
84    fOnlyWithSignal = kFALSE;
85    fThreshold = (Int_t)Conf()->GetNum("ChipTrackingTrgEngine0.fThreshold");
86
87    if ( Conf()->GetStr("ChipTrackingTrgEngine0.fThresholdType")=="absolute" ) fRelativeThreshold = kFALSE;
88    else if ( Conf()->GetStr("ChipTrackingTrgEngine0.fThresholdType")=="relative" ) fRelativeThreshold = kTRUE;
89    else Msg(EsafMsg::Panic) << "Wrong cfg value fThresholdType" << MsgDispatch;
90
91    ConfigFileParser *pConfig = Config::Get()->GetCF("Electronics","MacroCell");
92    if ( !pConfig->GetBool("MacroCell.fSaveAllChipGtuData") )
93        Msg(EsafMsg::Warning) << "MacroCell.fSaveAllChipGtuData is not enabled. Some infos are missing" << MsgDispatch;
94
95    FillRunPars();
96    Clear();
97}
98
99
100//______________________________________________________________________________
101ChipTrackingTrgEngine0::~ChipTrackingTrgEngine0() {
102    //
103    // Destructor
104    //
105
106    Clear();
107}
108
109
110//______________________________________________________________________________
111void ChipTrackingTrgEngine0::Clear() {
112    //
113    // Clear data containers
114    //
115
116    // clear track vectors
117    map<Int_t, vector<ChipTrackSegment> >::iterator it;
118    for(it=fTrackSegments.begin(); it!=fTrackSegments.end(); it++)
119        it->second.clear();
120
121    // clear map
122    fTrackSegments.clear();
123
124    for ( Int_t i(0); i<32; i++ ) nTracks[i] = 0;
125}
126
127//______________________________________________________________________________
128void ChipTrackingTrgEngine0::Simulate( MacroCellData* pData) {
129    //
130    // Simulate chip tracking trigger algorithm
131    //
132
133    fPatternCounter = 0;
134
135    // get number of channels in one chip
136    Int_t matsize = pData->Cell()->GetEC(0)->FrontEnd()->Channels();
137
138    // build a nearest neighbor matrix and its square
139    BoolMatrix A( matsize );
140    A.SetNeighborsMatrix();
141    BoolMatrix B( matsize );
142    B = A;
143    B *= A;
144
145    // get map of ChipGtuData elements
146    map<Int_t, map<Int_t,ChipGtuData*>* > &m1 = pData->fChipData2;
147    map<Int_t, map<Int_t,ChipGtuData*>* >::const_iterator it1;
148    // build 3 vectors of state vectors:
149    //   x vector is the list of state vectors v
150    //   y vector is the list of products A*v
151    //   z vector is the list of products A*(A*v)
152    vector<BoolVector> x;
153    vector<BoolVector> y;
154    vector<BoolVector> z;
155    vector<Int_t> gtu;
156
157    // iterate on chips in this macrocell
158    for(it1=m1.begin(); it1!=m1.end(); it1++) {
159        Int_t chip_id = it1->first;
160        map<Int_t,ChipGtuData*> &m2 = *(it1->second);
161        map<Int_t,ChipGtuData*>::const_iterator it2;
162        x.clear();
163        y.clear();
164        z.clear();
165        gtu.clear();
166
167        // iterate on GTU for the current chip
168        for(it2=m2.begin(); it2!=m2.end(); it2++) {
169            gtu.push_back(it2->first);
170        }
171
172        // if we do not have enough activity in this chip, skip to next one
173        if ( gtu.size() < (UInt_t)fMinTrackLength ) continue;
174
175        if (fOnlyWithSignal) {
176            Bool_t fIsSignal = kFALSE;
177            for(UInt_t i(0); i<gtu.size(); i++) {
178                fIsSignal += (Bool_t)(m2[gtu[i]]->GetTotalPureSignal());
179            }
180            if (!fIsSignal) continue;
181        }
182
183        Int_t mingtu = gtu[0];
184        Int_t maxgtu = gtu[0];
185        for ( UInt_t i(0); i<gtu.size(); i++) {
186            if ( gtu[i] < mingtu ) mingtu = gtu[i];
187            if ( gtu[i] > maxgtu ) maxgtu = gtu[i];
188        }
189        gtu.clear();
190
191        // compute the true threshold to use in trigger algorithm
192        // if fThresholdType=absolute, is simply the value specified in cfg file
193        Int_t fTrueThreshold;
194        if ( fRelativeThreshold ) {
195            Double_t gtu_length = Config::Get()->GetCF("Electronics","MacroCell")->GetNum("MacroCell.fGtuTimeLength");
196            // mean number of background hits in a pixel in a gtu
197            Double_t back = (m2[mingtu]->FrontEnd()->GetNightGlowRate()) * gtu_length;
198            //Msg(EsafMsg::Info) << "back " << back << " " << m2[mingtu]->FrontEnd()->GetNightGlowRate() << MsgDispatch;
199            fTrueThreshold = (Int_t)Ceil( back + fThreshold * Sqrt(back) );
200        } else {
201            fTrueThreshold = fThreshold;
202        }
203        //Msg(EsafMsg::Info) << "Chip " << chip_id << " thresh " << fTrueThreshold << MsgDispatch;
204        Int_t nPixelsOverThreshold(0);
205
206        for(Int_t i(0); i<(maxgtu-mingtu+1); i++) {
207            Int_t curr_gtu = mingtu + i;
208            gtu.push_back(curr_gtu);
209
210
211            // fill a boolean vector of pixels in this gtu
212            x.push_back( BoolVector(matsize) );
213            BoolVector& v1 = x.back();
214            v1.Zero();
215
216            if ( m2.count(curr_gtu) == 1 ) {
217                ChipGtuData &cgd = *(m2[curr_gtu]);
218                for(Int_t j=0; j<matsize; j++) {
219                    v1(j) = ( cgd.GetCounter(j) >= fTrueThreshold );    //true if j-th channel has at least fThreshold photoelectrons
220                    if (v1(j)) nPixelsOverThreshold++;
221                }
222            } else if (m2.count(curr_gtu)>1) {
223                Msg(EsafMsg::Warning) << "Duplicate ChipGtuData chip " << chip_id << " gtu " << curr_gtu << MsgDispatch;
224                ChipGtuData &cgd = *(m2[curr_gtu]);
225                for(Int_t j=0; j<matsize; j++) {
226                    v1(j) = ( cgd.GetCounter(j) >= fTrueThreshold );    //true if j-th channel has at least fThreshold photoelectrons
227                    if (v1(j)) nPixelsOverThreshold++;
228                }
229            }
230
231            y.push_back( BoolVector(matsize) );
232            BoolVector& v2 = y.back();
233            v2 = v1;
234            v2 *= A;
235
236            z.push_back( BoolVector(matsize) );
237            BoolVector& v3 = z.back();
238            v3 = v1;
239            v3 *= B;
240
241        }
242
243        nPixelsOverThreshold /= gtu.size();
244
245        // iterate on x and y to build tracks of length from minimum to maximum
246        // the algorythm accept at most one hole in the middle of any track
247        // (if this is required)
248
249        Int_t size = gtu.size();
250        Int_t last = size - fMinTrackLength + 1;
251        /*Msg(EsafMsg::Info) << "Number of valid GTU: "<< size << " last= "<< last << " Threshold: "
252                           << fThreshold << " Pixels over threshold: " << nPixelsOverThreshold << MsgDispatch;*/
253
254        //loop on valid gtus
255        for(Int_t i=0; i <= last ; i++) {
256
257            // loop on valid track length
258            for(Int_t k=(fMinTrackLength-1); k<=(fMaxTrackLength-1); k++) {
259
260                //Bool_t good_track = kFALSE;
261                Bool_t good_track2 = kFALSE;
262                Bool_t hole = kFALSE;
263
264                BoolVector trackVector(matsize);
265                trackVector.Zero();
266
267                // check if a track of length k can be found
268                if ( (i+k)>=size ) break;
269
270                // loop on hits of a possible track of length k
271                for(Int_t j=i; j < (i+k); j++) {
272
273                    // at the beginning there cannot be a hole
274                    if ( j==i ) {
275                        //good_track = (x[j])%(y[j+1]);
276                        trackVector = x[j];
277                        trackVector &= y[j+1];
278                        if ( !trackVector.OR() ) break;
279                    }
280
281                    // otherwise admit at most one hole
282                    else {
283                        //Bool_t a = x[j]%y[j+1];
284                        if ( x[j]%y[j+1] ) {
285                            //good_track &= a;
286                            trackVector &= y[j+1];
287                        }
288                        else if ( fAcceptHole && !hole && j<(i+k-1)) { //last point of tracks cannot be a hole
289                            //good_track &= x[j]%z[j+2];
290                            trackVector &= z[j+2];
291                            hole = kTRUE;
292                        }
293                        else {
294                            //good_track &= a;
295                            trackVector.Zero();
296                            break;
297                        }
298                    }
299                }
300
301                // good_track: original algorithm (boolean scalar product)
302                // good_track2: improved algorithm (bitwise AND + OR of the final vector)
303                good_track2 = trackVector.OR();
304
305                // if we succedeed to build a valid track, store it
306                if ( good_track2 ) {
307                    fTrackSegments[chip_id].push_back( ChipTrackSegment() );
308                    ChipTrackSegment &seg = fTrackSegments[chip_id].back();
309                    seg.SetGtuStart( gtu[i] );
310                    seg.SetGtuEnd( gtu[i+k] );
311                    seg.SetHasHole( hole );
312                    seg.SetTrackLength(k+1);
313                    seg.SetChipUid( chip_id );
314#ifdef DEBUG
315                    Msg(EsafMsg::Debug)<< "Track: len=" << seg.GetTrackLength() << " gtu:" << seg.GetGtuStart()
316                                       << "," << seg.GetGtuEnd() <<" chip_id: "<< chip_id << MsgDispatch;
317#endif /* DEBUG */
318                    // if track generate trigger, job is done
319                    if ( seg.GetTrackLength() >= fMinTriggerTrackLength ) {
320                        seg.SetTriggered( kTRUE );
321                        //SetTrigger( true );
322                        SetGtuTrigger( seg.GetGtuStart() );
323                        // check wich pattern of trigger is valid
324                        fPatternCounter = 0; // in trigger word first trigger in a microcell only
325                        for ( Int_t l=fMinTriggerTrackLength; l<=fMaxTrackLength; l++) {
326                            if (seg.GetTrackLength()==l) {
327                                SetTrigger( kTrigPattern[fPatternCounter] );
328                                nTracks[fPatternCounter]++;
329                            }
330                            fPatternCounter++;
331                        }
332                    }
333                }
334            }
335        }
336    }
337
338    // here try to attach tracks of different length
339    // in this version of the algorithm, we trigger if:
340    //    a track of length fMinTriggerTrackLnegth is found
341    //    2 tracks whose total sum is at least fMinTriggerTrackLnegth and the
342    //    end gtu of the first is connected to the first gtu of the second track
343    //    the second condition searches for tracks that are split in two ECs
344
345    for(it1=m1.begin(); it1!=m1.end(); it1++) {
346        Int_t chip_id = it1->first;
347        if ( !(fTrackSegments.count(chip_id))) continue;
348        map<Int_t,ChipGtuData*> &m2 = *(it1->second);
349        map<Int_t,ChipGtuData*>::const_iterator it2;
350        it2 = m2.begin();
351        Int_t uidSt = (*(it2->second)).FrontEnd()->UniqueChanId(0);
352        Int_t uidEnd = (*(it2->second)).FrontEnd()->UniqueChanId(matsize-1);
353        pair<Int_t,Int_t> rcSt = pData->Cell()->GetRowColUniqueId(uidSt);
354        pair<Int_t,Int_t> rcEnd = pData->Cell()->GetRowColUniqueId(uidEnd);
355        Int_t rmin(0),rmax(0),cmin(0),cmax(0);
356        if (rcSt.first > rcEnd.first ) {
357            rmax = rcSt.first;
358            rmin = rcEnd.first;
359        } else {
360            rmax = rcEnd.first;
361            rmin = rcSt.first;
362        }
363        if (rcSt.second > rcEnd.second ) {
364            cmax = rcSt.second;
365            cmin = rcEnd.second;
366        } else {
367            cmax = rcEnd.second;
368            cmin = rcSt.second;
369        }
370
371        // search the neighbors chips with at least a valid track segment
372        vector<Int_t> fChipNeighbors;
373        for(Int_t kk=0; kk<4; kk++) {
374            Int_t row0(0), col0(0), sr(0), sc(0);
375            switch(kk) {
376                case 0: {row0 = rmin; col0=cmin; sr=-1; sc=-1;}
377                case 1: {row0 = rmin; col0=cmax; sr=-1; sc=1;}
378                case 2: {row0 = rmax; col0=cmin; sr=1; sc=-1;}
379                case 3: {row0 = rmax; col0=cmax; sr=1; sc=1;}
380            }
381            for( Int_t r(0); r<2; r++ ) {
382                for( Int_t c(0); c<2; c++ ) {
383                    if ( r == 0 && c==0 ) continue;
384                    Int_t row = row0 + sr*r;
385                    Int_t col = col0 + sc*c;
386                    Int_t id = 0;
387                    Photomultiplier* pmt = GetEusoElectronics()->PmtId(pData->Cell()->GetUniqueIdRowCol(row,col));
388                    if (pmt) id = pmt->FrontEnd()->Id(); else continue;
389                    if (!id || id == chip_id || fTrackSegments.count(id) == 0) continue;
390                    Bool_t flag = kTRUE;
391                    for ( Int_t i(0); i<(Int_t)fChipNeighbors.size(); i++ ) {
392                        if ( fChipNeighbors[i] == id ) flag=kFALSE;
393                    }
394                    if (flag) fChipNeighbors.push_back(id);
395                }
396            }
397        }
398
399        // merge tracks segments between nearby chips
400        for(Int_t i(0); i<(Int_t)fTrackSegments[chip_id].size(); i++) {
401            ChipTrackSegment &cseg = fTrackSegments[chip_id][i];
402            //if ( cseg.GetTriggered()==kTRUE ) continue;
403            if ( cseg.GetTrackLength() > fMinTriggerTrackLength ) continue;    // only for my trigger studies (RP)
404            for( Int_t k(0); k<(Int_t)fChipNeighbors.size(); k++ ){
405                Int_t id = fChipNeighbors[k];
406                for( Int_t j(0); j<(Int_t)fTrackSegments[id].size(); j++ ) {
407                    ChipTrackSegment seg = fTrackSegments[id][j];
408                    Int_t totlen = seg.GetTrackLength() + cseg.GetTrackLength();
409                    if ( totlen < fMinTriggerTwoLength || totlen > fMaxTwoLength ) continue;
410                    if ( TMath::Abs(cseg.GetGtuEnd()-seg.GetGtuEnd()) < (seg.GetTrackLength()+1) ) {
411                        //cseg.SetTriggered(kTRUE);
412                        //SetTrigger( true );
413                        fPatternCounter = fMaxTrackLength - fMinTriggerTrackLength + 1; // skip patterns for a microcell only
414                        for(Int_t l=fMinTriggerTwoLength; l<fMaxTwoLength+1; l++) {
415                            if ( totlen == l ) {
416                                SetTrigger( kTrigPattern[fPatternCounter] );
417                                nTracks[fPatternCounter]++;
418                            }
419                            fPatternCounter++;
420                        }
421                        SetGtuTrigger( (seg.GetGtuStart() > cseg.GetGtuStart()) ? cseg.GetGtuStart() : seg.GetGtuStart() );
422                        if ( cseg.GetTriggered() == kFALSE ) {
423                            cseg.SetTriggered(kTRUE);
424                            fTrackSegments[chip_id][i] = cseg; }
425                        if ( seg.GetTriggered()==kFALSE ) {
426                            seg.SetTriggered(kTRUE);
427                            fTrackSegments[id][j] = seg;
428                        }
429                        //break;
430                    }
431                }
432                //if ( cseg.GetTriggered()==kTRUE ) break;
433            }
434        }
435        fChipNeighbors.clear();
436    }
437
438    Int_t lcounter = 0;
439    for ( Int_t i(0); i<(fMaxTrackLength-fMinTriggerTrackLength+1); i++ ) {
440        Msg(EsafMsg::Info) << "Found " << nTracks[lcounter] << " tracks of length " << fMinTriggerTrackLength + i << " in a single microcell" << MsgDispatch;
441        lcounter++;
442    }
443    for ( Int_t i(0); i<(fMaxTwoLength-fMinTriggerTwoLength+1); i++ ) {
444Msg(EsafMsg::Info) << "Found " << nTracks[lcounter] << " tracks of length " << fMinTriggerTwoLength + i << " in two adjacent microcells" << MsgDispatch;
445        lcounter++;
446    }
447
448    FillEEvent(pData->Cell()->Id());
449    Clear();
450}
451
452//_____________________________________________________________________________________
453void ChipTrackingTrgEngine0::FillEEvent( Int_t cell_id ) {
454    //
455    // Fill root event
456    //
457    map<Int_t, vector<ChipTrackSegment> >::iterator it;
458    for( it = fTrackSegments.begin(); it != fTrackSegments.end(); it++ ) {
459        Int_t chip_id = (*it).first;
460        for( UInt_t i(0); i<((*it).second).size(); i++ ) {
461            if ( EEvent::GetCurrent() ) {
462                EChipTrackTriggerDataAdder a( &(((*it).second)[i]), cell_id, chip_id );
463                EEvent::GetCurrent()->Fill(a);
464            }
465        }
466    }
467}
468
469//_____________________________________________________________________________________
470void ChipTrackingTrgEngine0::FillRunPars() {
471    //
472    // Fill run parameters
473    //
474
475    if (!(EEvent::GetCurrent() && EEvent::GetCurrent()->GetRunPars()) ) return;
476    ERunParameters *runpars = EEvent::GetCurrent()->GetRunPars();
477
478    if (runpars) {
479        EChipTriggParsFiller ctpfiller;
480        ctpfiller.SetName(GetName().c_str());
481        ctpfiller.SetId((ETriggerTypeIdentifier)GetTriggerId());
482        ctpfiller.SetRelativeThreshold(fRelativeThreshold);
483        ctpfiller.SetThreshold(fThreshold);
484        ctpfiller.SetMinTrackLength(fMinTrackLength);
485        ctpfiller.SetMaxTrackLength(fMaxTrackLength);
486        ctpfiller.SetMinTriggerTrackLength(fMinTriggerTrackLength);
487        ctpfiller.SetMinTriggerTwoLength(fMinTriggerTwoLength);
488        ctpfiller.SetMaxTwoLength(fMaxTwoLength);
489        ctpfiller.SetAcceptHole(fAcceptHole);
490        ctpfiller.SetOnlyWithSignal(fOnlyWithSignal);
491        runpars->Fill(ctpfiller);
492    }
493
494}
495
496//_____________________________________________________________________________________
497void ChipTrackingTrgEngine0::Dump( Int_t cell_id ) {
498    //
499    // Dump the number of tracks (triggered and total) for each chip
500    //
501    map<Int_t, vector<ChipTrackSegment> >::iterator it;
502    Msg(EsafMsg::Info) << "***DUMPING CELL " << cell_id << MsgDispatch;
503    if (fTrackSegments.size()==0) return;
504    for( it = fTrackSegments.begin(); it != fTrackSegments.end(); it++ ) {
505        Int_t chip_id = (*it).first;
506        Int_t tottracks = ((*it).second).size();
507        Int_t trgtracks(0);
508        Int_t gtu_min, gtu_max;
509        gtu_min = ((*it).second)[0].GetGtuStart();
510        gtu_max = ((*it).second)[0].GetGtuEnd();
511        for( UInt_t i(0); i<((*it).second).size(); i++ ) {
512            if ( ((*it).second)[i].GetTriggered() ) trgtracks++;
513            if ( ((*it).second)[i].GetGtuStart() < gtu_min ) gtu_min = ((*it).second)[i].GetGtuStart();
514            if ( ((*it).second)[i].GetGtuEnd() > gtu_max ) gtu_max = ((*it).second)[i].GetGtuEnd();
515        }
516        Msg(EsafMsg::Info) << "Chip " << chip_id << " Total = " << tottracks << " Triggered = " << trgtracks << MsgDispatch;
517        Msg(EsafMsg::Info) << "Gtu min = " << gtu_min << " Gtu max = " << gtu_max << MsgDispatch;
518    }
519}
Note: See TracBrowser for help on using the repository browser.