AlgFilterDigitListSR.cxx

Go to the documentation of this file.

// $Id: AlgFilterDigitListSR.cxx,v 1.6 2009/03/02 16:43:18 bckhouse Exp $
//
// AlgFilterDigitListSR.cxx
//
// Concrete Algorithm for filtering CandDigit's from CandDigitList
//
// Author:  R. Lee 2/28/2003
// AlgConfig parameters (with default values) are:
//
// "FilterStrategy" - Filtering strategy to use

#include <cassert>
#include <cmath>

#include "Algorithm/AlgConfig.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "CandData/CandHeader.h"
#include "CandData/CandRecord.h"
#include "CandDigit/CandDigit.h"
#include "CandDigit/CandDigitHandle.h"
#include "CandDigit/CandDigitList.h"
#include "CandDigit/CandDigitListHandle.h"
#include "Candidate/CandContext.h"
#include "Conventions/ReadoutType.h"
#include "FilterDigitSR/AlgFilterDigitListSR.h"
#include "MessageService/MsgService.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Plex/PlexHandle.h"
#include "Validity/VldContext.h"


ClassImp(AlgFilterDigitListSR)

//______________________________________________________________________
CVSID("$Id: AlgFilterDigitListSR.cxx,v 1.6 2009/03/02 16:43:18 bckhouse Exp $");

//______________________________________________________________________
AlgFilterDigitListSR::AlgFilterDigitListSR()
{
}

//______________________________________________________________________
AlgFilterDigitListSR::~AlgFilterDigitListSR()
{
}

static NavKey KeyFromDigitTime(const CandDigitHandle* digit) {
  Int_t itime = (Int_t)(digit->GetSubtractedTime()*1.e9);
  return itime;
}


//______________________________________________________________________
void AlgFilterDigitListSR::RunAlg(AlgConfig &ac, CandHandle &ch,
                                                        CandContext &cx)
{


   MSG("FilterDigitSR", Msg::kDebug)
                     << "Starting AlgFilterDigitListSR::RunAlg()" << endl;

// Check for input CandDigitListHandle in CandContext
   assert(cx.GetDataIn());

// Obtain CandContext's CandDigitList as input
   const CandDigitListHandle *cdlh =
             dynamic_cast<const CandDigitListHandle *> (cx.GetDataIn());
   assert(cdlh);   // Check for input CandDigitListHandle in CandContext

// Obtain filtering parameters from AlgConfig
   Int_t tmpi = 0;
   Double_t tmpd = 0.;
   Int_t filterstrategy = 0;
   Double_t dt_trig = 100.e-9;
   Double_t dt_readout = 150.e-9;
   Double_t dt_prereadout = 50.e-9;
   Int_t shieldmode = 0;
   if (ac.Get("FilterStrategy", tmpi)) filterstrategy = tmpi;
   if (ac.Get("TriggerWindow",tmpd)) dt_trig = (tmpd)*1.e-9;
   if (ac.Get("PostReadoutWindow",tmpd)) dt_readout = (tmpd)*1.e-9;
   if (ac.Get("PreReadoutWindow",tmpd)) dt_prereadout = (tmpd)*1.e-9;
   if (ac.Get("ShieldMode",tmpi)) shieldmode = tmpi;
   dt_trig = fabs(dt_trig);
   dt_readout = fabs(dt_readout);
   dt_prereadout = fabs(dt_prereadout);
   if (shieldmode<0 || shieldmode>2) shieldmode=0;

// Set CandRecord pointer in new CandDigitList to that of CandDigitList
   ch.SetCandRecord((cdlh->GetCandRecord())); //Sub-algorithms need this

// Copy member variables from orig CandDigitList to new CandDigitList
   CandDigitListHandle &cddlh = dynamic_cast<CandDigitListHandle &>(ch);
   cddlh.SetIsSparse(cdlh->GetIsSparse());
   cddlh.SetAbsTime(cdlh->GetAbsTime());

// Get VldContext
   const VldContext &vldc = *(cdlh->GetCandRecord()->GetVldContext());

   PlexHandle ph(vldc);

// Iterate over CandContext's CandDigitList to find CandDigitHandles.
   CandDigitHandleItr cdhiter(cdlh->GetDaughterIterator());
   CandDigitHandleKeyFunc *cdhKf = cdhiter.CreateKeyFunc();
   cdhKf->SetFun(KeyFromDigitTime);
   cdhiter.GetSet()->AdoptSortKeyFunc(cdhKf);
   cdhKf = 0;

   Float_t totadc = 0.;
   Float_t totadcfilter = 0.;

   Double_t endtime = 0.;
   TObject *tob;
   TObjArray digitlist;
   map<CandDigitHandle*,Bool_t> shielddigitmap;
   map<CandDigitHandle*,Bool_t>::iterator shielddigitmapiter;
   while ((tob = cdhiter())) {
     CandDigitHandle *cdh = dynamic_cast<CandDigitHandle *>(tob);
     if (cdh == 0) continue;   // Skip daughter if not a CandDigitHandle
     Bool_t isshielddigit = 0;
     if (cdh->GetPlexSEIdAltL().IsVetoShield()) isshielddigit=1;
     if (shieldmode==2 && isshielddigit) continue;
     if (isshielddigit) shielddigitmap[cdh] = 1;
     if (ph.GetReadoutType(cdh->GetChannelId()) != ReadoutType::kScintStrip) continue;
     if (!isshielddigit && cdh->GetSubtractedTime()<0.) continue; // pre digits not counted -- the 4/5 trigger calculation in the DAQ was done with digits GetSubtractedTime()>=0
     if (cdh->GetSubtractedTime()>endtime) endtime = cdh->GetSubtractedTime();
     totadc += cdh->GetCharge();
     digitlist.Add(cdh);
   }


   CandDigitHandle *gooddgt = 0;
   Float_t goodadc = 0.;
   Int_t timeplane[500]; // used to eliminate redundant calculations
   for (int i=0; i<500; i++) {
     timeplane[i] = -999999999;
   }
   for (int i=0; i<=digitlist.GetLast(); i++) {
     CandDigitHandle *dgt0 = dynamic_cast<CandDigitHandle*>(digitlist.At(i));
     shielddigitmapiter = shielddigitmap.find(dgt0);
     if (shielddigitmapiter!=shielddigitmap.end()) continue; // don't use shield digits for 4/5 calculation
     Float_t trigplane[500];
     for (int j=0; j<500; j++) trigplane[j] = 0.;
     Int_t iplane0 = dgt0->GetPlexSEIdAltL().GetPlane();
     if (iplane0<=0 ||iplane0>=500) continue;
     if (dgt0->GetSubtractedTime()<=timeplane[iplane0]) continue;
     timeplane[iplane0] = (Int_t)(dgt0->GetSubtractedTime()*1.e9);
     trigplane[iplane0] = dgt0->GetCharge();
     CandDigitHandle *dgt1 = 0;
     Int_t ipln0=iplane0;
     Int_t ipln1=iplane0;
     for (int j=i+1; j<=digitlist.GetLast() && (!dgt1 || dgt1->GetSubtractedTime()-dgt0->GetSubtractedTime()<=dt_trig); j++) {
       dgt1 = dynamic_cast<CandDigitHandle*>(digitlist.At(j));
       Int_t iplane1 = dgt1->GetPlexSEIdAltL().GetPlane();
       if (iplane1<=0 ||iplane1>=500) continue;
       trigplane[iplane1] += dgt1->GetCharge();
       if (iplane1<ipln0) ipln0 = iplane1;
       if (iplane1>ipln1) ipln1 = iplane1;
     }
     for (int iplane = ipln0; iplane<=ipln1-3; iplane++) {
       Int_t nhitplane = 0;
       Int_t mintrigplane = 9999;
       Int_t maxtrigplane = -1;
       Float_t thisadc = 0.;
       for (int dplane=0; dplane<5 && iplane+dplane<=ipln1; dplane++) {
         if (trigplane[iplane+dplane]>0.) {
           nhitplane++;
           thisadc += trigplane[iplane+dplane];
           if (iplane+dplane<mintrigplane) mintrigplane = iplane+dplane;
           if (iplane+dplane>maxtrigplane) maxtrigplane = iplane+dplane;
         }
       }
       if (nhitplane>=4 && iplane0>=mintrigplane && iplane0<=maxtrigplane && thisadc>goodadc) {
         gooddgt = dgt0;
         goodadc = thisadc;
       }
     }
   }

   if (gooddgt) {
     cddlh.SetAbsTime(gooddgt->GetTime());
     for (int i=0; i<=digitlist.GetLast(); i++) {
       CandDigitHandle *dgt = dynamic_cast<CandDigitHandle*>(digitlist.At(i));
       Bool_t isshielddigit = 0;
       shielddigitmapiter = shielddigitmap.find(dgt);
       if (shielddigitmapiter!=shielddigitmap.end()) isshielddigit = 1;
       if ((isshielddigit && shieldmode==0) || (gooddgt->GetSubtractedTime()-dgt->GetSubtractedTime()<=dt_prereadout && gooddgt->GetSubtractedTime()-dgt->GetSubtractedTime()>=0.) || (dgt->GetSubtractedTime()-gooddgt->GetSubtractedTime()<=dt_readout && dgt->GetSubtractedTime()-gooddgt->GetSubtractedTime()>=0.)) {
         CandDigitHandle cddh(*dgt);
         ch.AddDaughterLink(cddh, kFALSE);           // Don't check for dups
         totadcfilter += cddh.GetCharge();
         MSG("FilterDigitSR",Msg::kVerbose) << "adding digit, plane " << dgt->GetPlexSEIdAltL().GetPlane() << " relative time " << (dgt->GetSubtractedTime()-gooddgt->GetSubtractedTime())*1.e9 << " isshield " << isshielddigit << endl;
       }
     }
   }

   const CandHeader *candheader = cdlh->GetCandRecord()->GetCandHeader();
   assert(candheader);
   double fraction = 0;
   if (totadc != 0) fraction = totadcfilter/totadc;

   MSG("FilterDigitSR",Msg::kDebug) << "Run " << candheader->GetRun() << " Snarl " << candheader->GetSnarl() << " trigger = " << (gooddgt ? 1 : 0) << " adc = " << totadc << " fraction = " << fraction << " timelength = " << endtime*1.e9 << endl;


}

//______________________________________________________________________
void AlgFilterDigitListSR::Trace(const char * /* c */) const
{
}

---