Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

SignalChipTrackingTrgEngine - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: SignalChipTrackingTrgEngine.cc,v 1.4 2005/11/11 09:08:06 pesce Exp $
// R. Pesce created

//_____________________________________________________________________________
//  
//  Signal Chip Tracking Trigger Engine 
//  ===================================
//
//  Tracking trigger built in front end chip.
//  The algorithm searches for sequences of pixels (over the threshold)
//  consequtive and adjacent. If the track length is over a threshold,
//  trigger occurs.
//  The algorithm searches also between track segments of two nearby chips:
//  if two track segments are consequtive in time and the total length is
//  over a threshold (different from the previous one), trigger occurs
//
//  Only the pure signal is processed by this trigger engine
//
//  Parameters:
//
//  fThreshold: threshold on pixel hits in a gtu
//  fMinTrackLenght: min length of the valid tracks
//  fMaxTrackLength: max length of the valid tracks
//  fMinTriggerTrackLength: min track length to have a trigger (single chip)
//  fMinTriggerTwoLength: min total track length in 2 nearby chips
//  fMaxTwoLength: max track length allowed in 2 nearby chips
//
//  fAcceptHole [bool] : an hole in a single track segment is or not accepted
//
//  Special multiple configuration:
//  The alorithm counts the numebr of tracks of the various trigger lengths,
//  fill a specific trigger word and dump the numbers of the recognized tracks
//


#include "SignalChipTrackingTrgEngine.hh"
#include "MacroCellData.hh"
#include "MacroCell.hh"
#include "ElementaryCell.hh"
#include "FrontEndChip.hh"
#include "BoolAlgebra.hh"
#include "ChipTrackSegment.hh"
#include "EusoElectronics.hh"
#include "Photomultiplier.hh"
#include "EusoDetector.hh"
#include "TMath.h"
#include "ERunParameters.hh"
#include "EChipTrackTriggerDataAdder.hh"
#include "EChipTriggParsFiller.hh"
#include "Config.hh"

using namespace TMath;

// trigger pattern types
const Int_t kTrigPattern[32] = {BIT(0),BIT(1),BIT(2),BIT(3),BIT(4),BIT(5),BIT(6),BIT(7),
                                BIT(8),BIT(9),BIT(10),BIT(11),BIT(12),BIT(13),BIT(14),BIT(15),
                                BIT(16),BIT(17),BIT(18),BIT(19),BIT(20),BIT(21),BIT(22),BIT(23),
                                BIT(24),BIT(25),BIT(26),BIT(27),BIT(28),BIT(29),BIT(30),BIT(31)};


ClassImp(SignalChipTrackingTrgEngine)

//______________________________________________________________________________
 SignalChipTrackingTrgEngine::SignalChipTrackingTrgEngine() : 
TriggerEngine(string("SignalChipTrackingTrgEngine"), kSignalChipTrackingTrigger ) {
    //
    // Constructor
    //

    fMinTrackLength = (Int_t)Conf()->GetNum("SignalChipTrackingTrgEngine.fMinTrackLength");
    fMaxTrackLength = (Int_t)Conf()->GetNum("SignalChipTrackingTrgEngine.fMaxTrackLength");
    fMinTriggerTrackLength = (Int_t)Conf()->GetNum("SignalChipTrackingTrgEngine.fMinTriggerTrackLength");
    fMinTriggerTwoLength = (Int_t)Conf()->GetNum("SignalChipTrackingTrgEngine.fMinTriggerTwoLength");
    fMaxTwoLength = (Int_t)Conf()->GetNum("SignalChipTrackingTrgEngine.fMaxTwoLength");

    if ( Conf()->GetStr("SignalChipTrackingTrgEngine.fThresholdType")=="absolute" ) fRelativeThreshold = kFALSE;
    else if ( Conf()->GetStr("SignalChipTrackingTrgEngine.fThresholdType")=="relative" ) fRelativeThreshold = kTRUE;
    else Msg(EsafMsg::Panic) << "Wrong cfg value fThresholdType" << MsgDispatch;
    fAcceptHole = Conf()->GetBool("SignalChipTrackingTrgEngine.fAcceptHole");
    fOnlyWithSignal = kTRUE;
    fThreshold = (Int_t)Conf()->GetNum("SignalChipTrackingTrgEngine.fThreshold");

    // panic if we want to simulate >32 trigger patterns
    if ( ((fMaxTrackLength-fMinTriggerTrackLength+1)+(fMaxTwoLength-fMinTriggerTwoLength+1)) > 32  )
        Msg(EsafMsg::Panic) << "Request for > 32 patterns of trigger." << MsgDispatch;

    ConfigFileParser *pConfig = Config::Get()->GetCF("Electronics","MacroCell");
    if ( !pConfig->GetBool("MacroCell.fSaveAllChipGtuData") )
        Msg(EsafMsg::Warning) << "MacroCell.fSaveAllChipGtuData is not enabled. Some infos are missing" << MsgDispatch;

    FillRunPars();
    Clear();
}


//______________________________________________________________________________
 SignalChipTrackingTrgEngine::~SignalChipTrackingTrgEngine() {
    //
    // Destructor
    //
    
    Clear();
}


//______________________________________________________________________________
 void SignalChipTrackingTrgEngine::Clear() {
    //
    // Clear data containers
    //
    
    // clear track vectors
    map<Int_t, vector<ChipTrackSegment> >::iterator it;
    for(it=fTrackSegments.begin(); it!=fTrackSegments.end(); it++)
        it->second.clear();
    
    // clear map
    fTrackSegments.clear();

}

//______________________________________________________________________________
 void SignalChipTrackingTrgEngine::Simulate( MacroCellData* pData) {    
    //
    // Simulate chip tracking algorythm
    //

    fPatternCounter = 0;

#ifdef DEBUG    
    Msg(EsafMsg::Debug) << "=== CHIP TRACKING STARTED CELL " << pData->Cell()->Id() << MsgDispatch;
    Msg(EsafMsg::Debug) << "Rows: " << pData->Cell()->GetRows() << MsgDispatch;
    Msg(EsafMsg::Debug) << "Columns: " << pData->Cell()->GetColumns() << MsgDispatch;
    Msg(EsafMsg::Debug) << "Chips: " << pData->Cell()->Chips() << MsgDispatch;
#endif /* DEBUG */

    Int_t nTriggeredTracks = 0;

    // get number of channels in one chip
    Int_t matsize = pData->Cell()->GetEC(0)->FrontEnd()->Channels();

    // build a nearest neighbor matrix and its square
    BoolMatrix A( matsize );
    A.SetNeighborsMatrix();
    BoolMatrix B( matsize );
    B = A;
    B *= A;

    // get map of ChipGtuData elements
    map<Int_t, map<Int_t,ChipGtuData*>* > &m1 = pData->fChipData2;
    map<Int_t, map<Int_t,ChipGtuData*>* >::const_iterator it1;
#ifdef DEBUG
    if ( m1.size()!=0 ) Msg(EsafMsg::Info) << "MacroCell " << pData->Cell()->Id() << " has " << m1.size() << " active chips" << MsgDispatch;
#endif /* DEBUG */
    // build 3 vectors of state vectors: 
    //   x vector is the list of state vectors v
    //   y vector is the list of products A*v
    //   z vector is the list of products A*(A*v)
    vector<BoolVector> x;
    vector<BoolVector> y;
    vector<BoolVector> z;
    vector<Int_t> gtu;
    
    // iterate on chips in this macrocell
    for(it1=m1.begin(); it1!=m1.end(); it1++) {
        Int_t chip_id = it1->first;
#ifdef DEBUG 
        Msg(EsafMsg::Debug) << "Processing chip " << chip_id << MsgDispatch;
#endif /* DEBUG */

        map<Int_t,ChipGtuData*> &m2 = *(it1->second);
        map<Int_t,ChipGtuData*>::const_iterator it2;
#ifdef DEBUG 
        if ( m2.size()!=0 ) Msg(EsafMsg::Info) << "Chip " << chip_id << " has " << m2.size() << " chipgtudata" << MsgDispatch;
#endif /* DEBUG */

        x.clear();
        y.clear();
        z.clear();
        gtu.clear();

        // iterate on GTU for the current chip
        for(it2=m2.begin(); it2!=m2.end(); it2++) {
            gtu.push_back(it2->first);
        }

        // if we do not have enough activity in this chip, skip to next one
        if ( gtu.size() < (UInt_t)fMinTrackLength ) continue;
        
        if (fOnlyWithSignal) {
            Bool_t fIsSignal = kFALSE;
            for(UInt_t i(0); i<gtu.size(); i++) {
                fIsSignal += (Bool_t)(m2[gtu[i]]->GetTotalPureSignal());
            }
            if (!fIsSignal) continue;
        }
        
        Int_t mingtu = gtu[0];
        Int_t maxgtu = gtu[0];
        for ( UInt_t i(0); i<gtu.size(); i++) {
            if ( gtu[i] < mingtu ) mingtu = gtu[i];
            if ( gtu[i] > maxgtu ) maxgtu = gtu[i];
        }
        gtu.clear();

        // compute the true threshold to use in trigger algorithm
        // if fThresholdType=absolute, is simply the value specified in cfg file
        Int_t fTrueThreshold;
        if ( fRelativeThreshold ) {
            Double_t gtu_length = Config::Get()->GetCF("Electronics","MacroCell")->GetNum("MacroCell.fGtuTimeLength");
            // mean number of background hits in a pixel in a gtu
            Double_t back = (m2[mingtu]->FrontEnd()->GetNightGlowRate()) * gtu_length;
            //Msg(EsafMsg::Info) << "back " << back << " " << m2[mingtu]->FrontEnd()->GetNightGlowRate() << MsgDispatch; 
            fTrueThreshold = (Int_t)Ceil( back + fThreshold * Sqrt(back) );
        } else {
            fTrueThreshold = fThreshold;
        }
        //Msg(EsafMsg::Info) << "Chip " << chip_id << " thresh " << fTrueThreshold << MsgDispatch;
        
        for(Int_t i(0); i<(maxgtu-mingtu+1); i++) {
            Int_t curr_gtu = mingtu + i;
            gtu.push_back(curr_gtu);
            
            
            // fill a boolean vector of pixels in this gtu
            x.push_back( BoolVector(matsize) );
            BoolVector& v1 = x.back();
            v1.Zero();

            if ( m2.count(curr_gtu) == 1 ) {
                ChipGtuData &cgd = *(m2[curr_gtu]);
                for(Int_t j=0; j<matsize; j++)
                    // only pure signal is taken in account
                    v1(j) = ( cgd.GetPureSignal(j) >= fTrueThreshold );    //true if j-th channel has at least fThreshold photoelectrons
                                                                    //of pure signal
            } else if (m2.count(curr_gtu)>1) {
                Msg(EsafMsg::Warning) << "Duplicate ChipGtuData chip " << chip_id << " gtu " << curr_gtu << MsgDispatch;
                ChipGtuData &cgd = *(m2[curr_gtu]);
                for(Int_t j=0; j<matsize; j++)
                    // only pure signal is taken in account
                    v1(j) = ( cgd.GetPureSignal(j) >= fTrueThreshold );    //true if j-th channel has at least fThreshold photoelectrons
                                                                    //of pure signal
               }
            
            y.push_back( BoolVector(matsize) );
            BoolVector& v2 = y.back();
            v2 = v1;
            v2 *= A;
                        
            z.push_back( BoolVector(matsize) );
            BoolVector& v3 = z.back();
            v3 = v1;
            v3 *= B;
                        
        }


        // iterate on x and y to build tracks of length from minimum to maximum
        // the algorythm accept at most one hole in the middle of any track
        // (if this is required)

        Int_t size = gtu.size();
        Int_t last = size - fMinTrackLength + 1;
#ifdef DEBUG
        Msg(EsafMsg::Debug) << "Number of valid GTU: "<< size << "last= "<<last<< MsgDispatch;
#endif /* DEBUG */

        // loop on valid gtus
        for(Int_t i=0; i <= last ; i++) {

            // loop on valid track length
            for(Int_t k=(fMinTrackLength-1); k<=(fMaxTrackLength-1); k++) {

                //Bool_t good_track = kFALSE;
                Bool_t good_track2 = kFALSE;
                Bool_t hole = kFALSE;

                BoolVector trackVector(matsize);
                trackVector.Zero();
  
                // check if a track of length k can be found
                if ( (i+k)>=size ) break;
                                
                // loop on hits of a possible track of length k
                for(Int_t j=i; j < (i+k); j++) {
                    
                    // at the beginning there cannot be a hole
                    if ( j==i ) {
                        //good_track = (x[j])%(y[j+1]);
                        trackVector = x[j];
                        trackVector &= y[j+1];
                        if ( !trackVector.OR() ) break;
                    }

                    // otherwise admit at most one hole 
                    else {
                        //Bool_t a = x[j]%y[j+1];
                        if ( x[j]%y[j+1] ) {
                            //good_track &= a;
                            trackVector &= y[j+1];
                        }
                        else if ( fAcceptHole && !hole && j<(i+k-1)) {  //last point of tracks cannot be a hole
                            //good_track &= x[j]%z[j+2];
                            trackVector &= z[j+2];
                            hole = kTRUE;
                        }
                        else {
                            //good_track &= a;
                            trackVector.Zero();
                            break;
                        }
                    }
                }
               
                // good_track: original algorithm (boolean scalar product)
                // dood_track2: improved algorithm (bitwise AND + OR of the final vector)
                good_track2 = trackVector.OR();              

                
                // if we succedeed to build a valid track, store it
                if ( good_track2 ) {
                    fTrackSegments[chip_id].push_back( ChipTrackSegment() );
                    ChipTrackSegment &seg = fTrackSegments[chip_id].back();
                    seg.SetGtuStart( gtu[i] );
                    seg.SetGtuEnd( gtu[i+k] );
                    seg.SetHasHole( hole );
                    seg.SetTrackLength(k+1);
                    seg.SetChipUid( chip_id );
#ifdef DEBUG
                    Msg(EsafMsg::Debug)<< "Track: len=" << seg.GetTrackLength() << " gtu:" << seg.GetGtuStart()
                                       << "," << seg.GetGtuEnd() <<" chip_id: "<< chip_id << MsgDispatch;
#endif /* DEBUG */
                    // if track generate trigger, job is done
                    if ( seg.GetTrackLength() >= fMinTriggerTrackLength ) {
                        seg.SetTriggered( kTRUE );
                        //SetTrigger( true );
                        SetGtuTrigger( seg.GetGtuStart() );
                        // check wich pattern of trigger is valid
                        fPatternCounter = 0; // in trigger word first trigger in a microcell only
                        for ( Int_t l=fMinTriggerTrackLength; l<=fMaxTrackLength; l++) {
                            if (seg.GetTrackLength()==l) SetTrigger( kTrigPattern[fPatternCounter] );
                            fPatternCounter++;
                        }
                        nTriggeredTracks++;
                    }
                }
            }
        }
    }

    // here try to attach tracks of different length
    // in this version of the algorythm, we trigger if:
    //    a track of length fMinTriggerTrackLnegth is found
    //    2 tracks whose total sum is at least fMinTriggerTwoLength and the
    //    end gtu of the first is connected to the first gtu of the second track
    //    the second condition searches for tracks that are split in two ECs

    for(it1=m1.begin(); it1!=m1.end(); it1++) {
        Int_t chip_id = it1->first;
        if ( !(fTrackSegments.count(chip_id))) continue;
        map<Int_t,ChipGtuData*> &m2 = *(it1->second);
        map<Int_t,ChipGtuData*>::const_iterator it2;
        it2 = m2.begin();
        Int_t uidSt = (*(it2->second)).FrontEnd()->UniqueChanId(0);
        Int_t uidEnd = (*(it2->second)).FrontEnd()->UniqueChanId(matsize-1);
        pair<Int_t,Int_t> rcSt = pData->Cell()->GetRowColUniqueId(uidSt);
        pair<Int_t,Int_t> rcEnd = pData->Cell()->GetRowColUniqueId(uidEnd);
        Int_t rmin(0),rmax(0),cmin(0),cmax(0);
        if (rcSt.first > rcEnd.first ) {
            rmax = rcSt.first;
            rmin = rcEnd.first;
        } else {
            rmax = rcEnd.first;
            rmin = rcSt.first;
        }
        if (rcSt.second > rcEnd.second ) {
            cmax = rcSt.second;
            cmin = rcEnd.second;
        } else {
            cmax = rcEnd.second;
            cmin = rcSt.second;
        }

        // search the neighbors chips with at least a valid track segment
        vector<Int_t> fChipNeighbors;
        for(Int_t kk=0; kk<4; kk++) {
            Int_t row0(0), col0(0), sr(0), sc(0);
            switch(kk) {
                case 0: {row0 = rmin; col0=cmin; sr=-1; sc=-1;}
                case 1: {row0 = rmin; col0=cmax; sr=-1; sc=1;}
                case 2: {row0 = rmax; col0=cmin; sr=1; sc=-1;}
                case 3: {row0 = rmax; col0=cmax; sr=1; sc=1;}
            }
            for( Int_t r(0); r<2; r++ ) {
                for( Int_t c(0); c<2; c++ ) {
                    if ( r == 0 && c==0 ) continue;
                    Int_t row = row0 + sr*r;
                    Int_t col = col0 + sc*c;
                    Int_t id = 0;
                    Photomultiplier* pmt = GetEusoDetector()->GetEusoElectronics()->PmtId(pData->Cell()->GetUniqueIdRowCol(row,col));
                    if (pmt) id = pmt->FrontEnd()->Id(); else continue;
                    if (!id || id == chip_id || fTrackSegments.count(id) == 0) continue;
                    Bool_t flag = kTRUE;
                    for ( Int_t i(0); i<(Int_t)fChipNeighbors.size(); i++ ) {
                        if ( fChipNeighbors[i] == id ) flag=kFALSE;
                    }                    
                    if (flag) fChipNeighbors.push_back(id);
                }
            }
        }

        // merge tracks segments between nearby chips
        for(Int_t i(0); i<(Int_t)fTrackSegments[chip_id].size(); i++) {
            ChipTrackSegment &cseg = fTrackSegments[chip_id][i];
            //if ( cseg.GetTriggered()==kTRUE ) continue;
            if ( cseg.GetTrackLength() > fMinTriggerTrackLength ) continue;    // only for my trigger studies (RP)
            for( Int_t k(0); k<(Int_t)fChipNeighbors.size(); k++ ){
                Int_t id = fChipNeighbors[k];
                for( Int_t j(0); j<(Int_t)fTrackSegments[id].size(); j++ ) {
                    ChipTrackSegment seg = fTrackSegments[id][j];
                    Int_t totlen = seg.GetTrackLength() + cseg.GetTrackLength();
                    if ( totlen < fMinTriggerTwoLength || totlen > fMaxTwoLength ) continue;
                    if ( TMath::Abs(cseg.GetGtuEnd()-seg.GetGtuEnd()) < (seg.GetTrackLength()+1) ) {
                        //cseg.SetTriggered(kTRUE);
                        //SetTrigger( true );
                        // check which pattern of trigger and set the trigger word
                        fPatternCounter = fMaxTrackLength - fMinTriggerTrackLength + 1; // skip patterns for a microcell only
                        for(Int_t l=fMinTriggerTwoLength; l<fMaxTwoLength; l++) {
                            if ( totlen == l ) SetTrigger( kTrigPattern[fPatternCounter] );
                            fPatternCounter++;
                        }
                        SetGtuTrigger( (seg.GetGtuStart() > cseg.GetGtuStart()) ? cseg.GetGtuStart() : seg.GetGtuStart() );
                        if ( cseg.GetTriggered() == kFALSE ) { 
                            cseg.SetTriggered(kTRUE); 
                            nTriggeredTracks++;
                            fTrackSegments[chip_id][i] = cseg; }
                        if ( seg.GetTriggered()==kFALSE ) {
                            seg.SetTriggered(kTRUE);
                            nTriggeredTracks++;
                            fTrackSegments[id][j] = seg;
                        }
                        //break;
                    }
                }
                //if ( cseg.GetTriggered()==kTRUE ) break;
            }
        }
        fChipNeighbors.clear();
    }
    //if ( nTriggeredTracks )
        //Msg(EsafMsg::Info) << nTriggeredTracks << " tracks have triggered in cell " << pData->Cell()->Id() << MsgDispatch;

    //FillEEvent(pData->Cell()->Id());
    //Dump(pData->Cell()->Id());
    Clear();
}

//_____________________________________________________________________________________
 void SignalChipTrackingTrgEngine::FillEEvent( Int_t cell_id ) {
    //
    // Fill root event
    //
    map<Int_t, vector<ChipTrackSegment> >::iterator it;
    for( it = fTrackSegments.begin(); it != fTrackSegments.end(); it++ ) {
        Int_t chip_id = (*it).first;
        for( UInt_t i(0); i<((*it).second).size(); i++ ) {
            if ( EEvent::GetCurrent() ) {
                EChipTrackTriggerDataAdder a( &(((*it).second)[i]), cell_id, chip_id );
                EEvent::GetCurrent()->Fill(a);
            }
        }
    }
}

//_____________________________________________________________________________________
 void SignalChipTrackingTrgEngine::FillRunPars() {
    //
    // Fill run parameters
    //

    if (!(EEvent::GetCurrent() && EEvent::GetCurrent()->GetRunPars()) ) return;
    ERunParameters *runpars = EEvent::GetCurrent()->GetRunPars();

    if (runpars) {
        EChipTriggParsFiller ctpfiller;
        ctpfiller.SetName(GetName().c_str());
        ctpfiller.SetId((ETriggerTypeIdentifier)GetTriggerId());
        ctpfiller.SetRelativeThreshold(fRelativeThreshold);
        ctpfiller.SetThreshold(fThreshold);
        ctpfiller.SetMinTrackLength(fMinTrackLength);
        ctpfiller.SetMaxTrackLength(fMaxTrackLength);
        ctpfiller.SetMinTriggerTrackLength(fMinTriggerTrackLength);
        ctpfiller.SetMinTriggerTwoLength(fMinTriggerTwoLength);
        ctpfiller.SetMaxTwoLength(fMaxTwoLength);
        ctpfiller.SetAcceptHole(fAcceptHole);
        ctpfiller.SetOnlyWithSignal(fOnlyWithSignal);
        runpars->Fill(ctpfiller);
    }

}

//_____________________________________________________________________________________
 void SignalChipTrackingTrgEngine::Dump( Int_t cell_id ) {
    //
    // Dump the number of tracks (triggered and total) for each chip
    //
    map<Int_t, vector<ChipTrackSegment> >::iterator it;
    Msg(EsafMsg::Info) << "***DUMPING CELL " << cell_id << MsgDispatch;
    if (fTrackSegments.size()==0) return;
    for( it = fTrackSegments.begin(); it != fTrackSegments.end(); it++ ) {
        Int_t chip_id = (*it).first;
        Int_t tottracks = ((*it).second).size();
        Int_t trgtracks(0);
        Int_t gtu_min, gtu_max;
        gtu_min = ((*it).second)[0].GetGtuStart(); 
        gtu_max = ((*it).second)[0].GetGtuEnd();
        for( UInt_t i(0); i<((*it).second).size(); i++ ) {
            if ( ((*it).second)[i].GetTriggered() ) trgtracks++;
            if ( ((*it).second)[i].GetGtuStart() < gtu_min ) gtu_min = ((*it).second)[i].GetGtuStart();
            if ( ((*it).second)[i].GetGtuEnd() > gtu_max ) gtu_max = ((*it).second)[i].GetGtuEnd();
        }
        Msg(EsafMsg::Info) << "Chip " << chip_id << " Total = " << tottracks << " Triggered = " << trgtracks << MsgDispatch;
        Msg(EsafMsg::Info) << "Gtu min = " << gtu_min << " Gtu max = " << gtu_max << MsgDispatch;
    }
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005