SimDetector.cxx

Go to the documentation of this file.

#include "SimDetector.h"
#include "SimPmt.h"
#include "SimPmtMaker.h"

#include "SimVaElectronics.h"
#include "SimElectronics.h"
#include "SimElecMaker.h"

#include "SimDaqTrigger.h"
#include "SimDaqTrigMaker.h"

#include "Digitization/DigiScintHit.h"
#include "Digitization/DigiPE.h"

#include "MessageService/MsgService.h"

SimDetector* gDetector = NULL;

CVSID("$Id: SimDetector.cxx,v 1.58 2008/11/18 17:31:02 rhatcher Exp $");
ClassImp(SimDetector)

SimDetector::SimDetector( VldContext context, TRandom* random ) :
  fContext(context),
  fPlex(context),
  fRandom(random),
  fFarElectronics(NULL),
  fNearElectronics(NULL)
{
  if(fRandom ==NULL) {
    fRandom = gRandom;
  }

  // Hardwire the configuration defaults to the safest.
  fClearPmtList = 0;
  fM64Model = "SimPmtM64Full";
  fM16Model = "SimPmtUTM16";
  fFarElecModel  = "SimVaElectronics";
  fNearElecModel = "SimQieElectronics";

  // Electronics.
  fFarElectronics  = SimElecMaker::Create(fFarElecModel,context,fRandom);
  fNearElectronics = SimElecMaker::Create(fNearElecModel,context,fRandom);
  
  fDaqTriggerMask = 0xFF;
  for(int i=0;i<kSimMaxTriggers;i++)
    fDaqTriggers[i] = 0;  
  
  MSG("DetSim",Msg::kDebug) << "SimDetector::Constructor. Context: " << fContext.AsString() << endl;

  gDetector = this;
}

SimDetector::~SimDetector( void )
{
  if(fFarElectronics) delete fFarElectronics;
  if(fNearElectronics) delete fNearElectronics;

  // Wipe digit list.
  fDigits.clear();

  // Wipe pmt list.
  SimPmtList::iterator it;
  for(it = fPmts.begin(); it!= fPmts.end(); it++)  {
    if(it->second) delete it->second;
    fPmts.erase(it);
  }
  fPmts.clear();
  fActivePmts.clear();
}

Registry& SimDetector::DefaultConfig( void )
{
   // Supply default configuration.
  static Registry r; // Default configuration for module
    
  // Set name of config
  std::string name = "SimDetector.config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();

  r.Set("clearPmtList",0);

  // PMT model
  r.Set("M64Model","SimPmtM64Oxford");
  r.Set("M16Model","SimPmtM16UTTimed");
  r.Set("pmtM16DefaultGain",60);  // Used by decalibrator
  r.Set("pmtM64DefaultGain",100); // ditto
  r.Set("pmtDoOpticalCrosstalk",1);
  r.Set("pmtDoChargeCrosstalk",1);
  r.Set("pmtDoNonlinearity",0);
  r.Set("pmtDoDarkNoise",1);
  r.Set("pmtApplyGainDrift",0); // 0 = Trust PMT gains.
  r.Set("pmtScaleOpticalCrosstalk",1.0);
  r.Set("pmtScaleChargeCrosstalk",1.0);
  //Added by G. Pawloski to scale Adjacent and Diagonal separately
  r.Set("pmtScaleAdjacentChargeCrosstalk",1.0);
  r.Set("pmtScaleDiagonalChargeCrosstalk",1.0);
  r.Set("pmtScaleAdjacentOpticalCrosstalk",1.0);
  r.Set("pmtScaleDiagonalOpticalCrosstalk",1.0);
  r.Set("pmtHamamatsuPixelNumbering",0);
  r.Set("pmtChargeSlop",1.6);
  r.Set("pmtDoChargeSmear",1);
  r.Set("pmtM64DynodeSkipRate",0.00);
  r.Set("pmtM64NLThreshold",3.0);
  r.Set("pmtM64NonlinearityScale",1.0);
  r.Set("pmtRebuildGainMap",0);

  // Electronics.
  r.Set("farElecModel", "SimVaTimedElectronics");
  r.Set("nearElecModel","SimQieElectronics");
  r.Set("doLookupNonlinearity",1);

  // FD front end
  r.Set("vaGain", 0.375/Munits::fC );
  r.Set("vaPedWidth", 2.8);
  r.Set("vaTimeSmearingMode", SimTimeSmearingMode::kBreitWigner);
  r.Set("vaTimingWidth", 0.);
  r.Set("vaSparsifyThresh", 17);
  r.Set("vaDynodeThresh", 57);
  r.Set("vaNonlinearity",1);
  r.Set("vaLinearityP0",24.);
  r.Set("vaLinearityP1",16091.3);
  r.Set("vaLinearityP2",4.63);
  r.Set("vaLinearityP3",13267.);
  r.Set("vaLinearityP4",7085.12);
    
  // FD back end
  r.Set("varcTriggerMode", SimVarcTriggerMode::k2of36); // 2 out of 36.
  r.Set("varcTriggerWindow", 400.*Munits::ns);
  r.Set("vaChipRandomDeadRate",700.*Munits::hertz);
  r.Set("vaChipDeadTime",5*Munits::microsecond);
  r.Set("vaShapingTime",500.*Munits::ns);

  // ND front end.
  r.Set("qieDoLookukp",1); // Turn on lookup table.
  r.Set("qieDacCharge",1.4*Munits::fC); // Charge for one DAC count. Charge is in Munits (Coulombs)
  r.Set("qieDaqPerRangeZeroAac",1.6); // DAC counts per ADC in range 0
  r.Set("qieFlipLowAdc",30.);         // Point at which the range flips    
  r.Set("qieFlipHighAdc",180.);       // Point at which the range flips
  r.Set("qiePedestalWidthAdc",1.5);   // RMS of pedestal in ADC counts
  r.Set("qieAdcRms",0.0);             // RMS of a large charge in ADC counts
  r.Set("qiePedestalDac",50);         // Pedestal in DAC counts.
  r.Set("qieSparsifyThresh",70);      // Sparsification threshold (should be larger than pedestal!)

  // Daq trigger
  r.Set("applyTrigger",3);
  r.Set("daqTrigTimingMethod",SimDaqTrigger::kGap);
  r.Set("daqTrigTimeGap",156.0*Munits::ns); 
  r.Set("daqTrigTimeWindow",117.0*Munits::ns); 
  r.Set("daqTrigReadoutWindowStart",5500.*Munits::ns); // Pretrigger window, ns
  r.Set("daqTrigReadoutWindowStop",187*Munits::ns); // Posttrigger window, ns
  r.Set("daqTriggerMask",0xFF);

  // Bit  2 (hex 0x04)
  // Standard 4 out of 5 trigger
  Registry trigger2;
  trigger2.SetName("trigger2");
  trigger2.Set("type","PlaneTrigger");
  trigger2.Set("thresh",0.0);
  trigger2.Set("N",4);
  trigger2.Set("M",5);     
  r.Set("trigger2",trigger2);

  // Bit 5 (hex 0x20)
  // Standard energy trigger:
  // In a 4-plane window, demand
  // >=1500 ADCs, planes hit >= 2, N_hits >= 6
  Registry trigger5;
  trigger5.SetName("trigger5");
  trigger5.Set("type","EnergyTrigger");
  trigger5.Set("thresh",1500.);
  trigger5.Set("Nhits",6);
  trigger5.Set("N",2);
  trigger5.Set("M",4);     
  r.Set("trigger5",trigger5);

  // Add a whole bunch more triggers to be filled by the user.
  // This adds empty registries that can easily be modified.
  for(int i=0;i<kSimMaxTriggers;i++) {
    Registry triggerN;
    if(!(r.Get(Form("trigger%d",i),triggerN))) {
      triggerN.SetName(Form("trigger%d",i));
      r.Set(Form("trigger%d",i),triggerN);
    }
  }
  
  return r;
}

void  SimDetector::Config( Registry& config )
{
  // Modify the configuration.
  const char* s;

  //-- Operational concerns.
  config.Get("clearPmtList",fClearPmtList);

  //-- Pmt models.
  config.Get("M64Model",s);
  if(SimPmtMaker::HasRegisteredClass(s)) fM64Model = s;
  else 
    MSG("DetSim",Msg::kError) << "M64 Pmt Model " << s 
                              << " is not known. Reverting to default." << endl;
  
  config.Get("M16Model",s); fM16Model = s;
  if(SimPmtMaker::HasRegisteredClass(s)) fM16Model = s;
  else 
    MSG("DetSim",Msg::kError) << "M16 Pmt Model " << s 
                              << " is not known. Reverting to default." << endl;
  
  // Create dummy versions of the PMT to configure the model (with static functions)
  SimPmt* m16 = SimPmtMaker::Create(fM16Model,PlexPixelSpotId(),
                                    fContext,fRandom);
  m16->Config(config);
  delete m16;
  
  SimPmt* m64 = SimPmtMaker::Create(fM64Model,PlexPixelSpotId(),
                                    fContext,fRandom);
  m64->Config(config);
  delete m64;

  //-- Recreate and configure the electronics models.
  config.Get("farElecModel",s);
  if(SimElecMaker::HasRegisteredClass(s)) fFarElecModel = s;
  else 
    MSG("DetSim",Msg::kError) << "Far Detector Electronics Model " << s 
                              << " is not known. Reverting to default."  << endl;
  config.Get("nearElecModel",s);
  if(SimElecMaker::HasRegisteredClass(s)) fNearElecModel = s;
  else 
    MSG("DetSim",Msg::kError) << "Near Detector Electronics Model " << s 
                              << " is not known. Reverting to default."  << endl;

  // Gentlemen, we can rebuild him..
  if(fFarElectronics)  delete fFarElectronics;
  if(fNearElectronics) delete fNearElectronics;
  fFarElectronics  = SimElecMaker::Create(fFarElecModel, fContext,fRandom);
  fNearElectronics = SimElecMaker::Create(fNearElecModel,fContext,fRandom);
  // And boot them.
  fFarElectronics ->Config(config);
  fNearElectronics->Config(config);

  //-- Triggers.
  config.Get("applyTrigger",             fApplyTrigger);
  config.Get("daqTrigTimingMethod",      fDaqTrigTimingMethod);
  config.Get("daqTrigTimeGap",           fDaqTrigTimeGap);
  config.Get("daqTrigTimeWindow",        fDaqTrigTimeWindow);
  config.Get("daqTrigReadoutWindowStart",fDaqTrigReadoutWindowStart);
  config.Get("daqTrigReadoutWindowStop", fDaqTrigReadoutWindowStop);
  config.Get("daqTriggerMask",           fDaqTriggerMask);

  // Rebuild trigger list.
  Registry trigConfig;
  for(int itrig = 0; itrig<kSimMaxTriggers; itrig++) {
    if(fDaqTriggers[itrig]) delete fDaqTriggers[itrig]; // Delete old ones.
    fDaqTriggers[itrig] = 0;
    
    if( config.Get(Form("trigger%d",itrig),trigConfig) ) {
      const char* triggerType = 0;
      if( trigConfig.Get("type",triggerType) ) {
        fDaqTriggers[itrig] = SimDaqTrigMaker::Create(triggerType);
        if(fDaqTriggers[itrig]) {
          fDaqTriggers[itrig]->Configure(trigConfig);     
        } else {
          MSG("DetSim",Msg::kError) << "Trigger " << itrig << ": type=" << triggerType << " not known!" << endl;
        }
      } 
    }
  }
  
}


void SimDetector::Reset( VldContext context )
{
  MSG("DetSim",Msg::kDebug) << "SimDetector::Reset. Context: " << fContext.AsString() << endl;

  // Wipe digit list.
  fDigits.clear();
  fSnarls.clear();

  // Reset the context.
  fContext = context;
  fPlex = PlexHandle(context); 


  // Delete all PMTs in list: this means that 'duds' won't eat random numbers.
  if(fClearPmtList) {
    // The nasty way: delete every PMT every event.
    // This is slower, but safer for event-reproducability and heisenbug close encounters
    SimPmtList::iterator it;
    for(it = fPmts.begin(); it!= fPmts.end(); it++)  {
      if(it->second) delete it->second;
      fPmts.erase(it);
    }
    fPmts.clear();
  } else {
    // Reset all pmts in list.
    // Much more efficient and aesthetically pleasing.
    SimPmtList::iterator it;
    for(it = fPmts.begin(); it!= fPmts.end(); it++)  {
      SimPmt* pmt = it->second;
      pmt->Reset(context);
    }
  }
   
  fActivePmts.clear(); // Wipe list of active pmts.

  if(SimPixelTimeBucket::GetGlobalPixelBuckets()>0) {
    MSG("DetSim",Msg::kError) << "Error: Memory leak detected: total SimPixelTimeBuckets after reset is " 
                              << SimPixelTimeBucket::GetGlobalPixelBuckets() << endl;
  }

  // Reset the Electronics.
  if(fNearElectronics) fNearElectronics->Reset(fContext);
  if(fFarElectronics)  fFarElectronics->Reset(fContext);
  
  // Reset temporary lists.
  fDigiPeStorage.clear();
}

void SimDetector::Clean()
{
  // Wipe digit list.
  fDigits.clear();
  
  // Wipe pmt list.
  SimPmtList::iterator it;
  for(it = fPmts.begin(); it!= fPmts.end(); it++)  {
    if(it->second) delete it->second;
    fPmts.erase(it);
  }
  fPmts.clear();
}


void SimDetector::AddDigiPEList( const TObjArray* peListObj ) 
{
  
  fPeList = peListObj;
  if(!fPeList) MSG("DetSim",Msg::kError) << "Null PE list input into SimDetector.. will probably crash." << endl;
}





void SimDetector::Simulate( void )
{
  DistributePEs();
  SimulatePmts();
  SimulateElectronics();
  SimulateTrigger();
  CompileStats();
}

int SimDetector::GetNumberOfDigits(void)
{
  return fDigits.size();
}


SimDigit* SimDetector::GetSimDigit(UInt_t ihit)
{
  if(ihit<fDigits.size()) 
    return &(fDigits[ihit]);
  else return NULL;
}

void SimDetector::Print(Option_t* option) const
{
  std::string opt = option;
  Bool_t verbose = (opt.find("v") != std::string::npos);;
  Bool_t showpmts = (opt.find("p") != std::string::npos);;
  Bool_t showdigits = (opt.find("d") != std::string::npos);;
  Bool_t showsnarldigits = (opt.find("s") != std::string::npos);;
  
  if(verbose) {
    showpmts = showdigits = showsnarldigits = true;
  }

  printf("SimDetector:");
  printf("Current context: %s\n",fContext.AsString());
  printf("---------Config------------------------------\n");
  printf("PMT sim:       M16Model          %s\n",fM16Model.c_str());
  printf("               M64Model          %s\n",fM64Model.c_str());
  SimPmt::PrintConfig();
  if(fFarElectronics)  fFarElectronics->Print();
  if(fNearElectronics) fNearElectronics->Print();
  printf("Triggers:      fDaqTriggerMask            0x%x\n",fDaqTriggerMask);
  printf("               fDaqTrigTimingMethod       %s\n",
         (fDaqTrigTimingMethod==SimDaqTrigger::kGap)?"Gap":"Window");
  printf("               fDaqTrigTimeGap            %f ns\n",fDaqTrigTimeGap/Munits::ns);
  printf("               fDaqTrigTimeWindow         %f ns\n",fDaqTrigTimeWindow/Munits::ns);
  printf("               fDaqTrigReadoutWindowStart %f ns\n",fDaqTrigReadoutWindowStart/Munits::ns);
  printf("               fDaqTrigReadoutWindowStop  %f ns\n",fDaqTrigReadoutWindowStop/Munits::ns);

  for(int itrig=0;itrig<kSimMaxTriggers;itrig++) {
    if(fDaqTriggers[itrig]) {
     printf("   -- Trigger %2d: --\n",itrig);
     fDaqTriggers[itrig]->Print(option);
    }
  }
  printf("---------PMTS------------------------------------\n");
  printf("    Total PMTS in memory: %d  Active PMTs: %d\n", (int)fPmts.size(),(int)fActivePmts.size());
  SimPmtList::iterator it;
  if(showpmts)
    for(it = fActivePmts.begin(); it!= fActivePmts.end(); it++)  {
      SimPmt* pmt = it->second;
      if(pmt){
        pmt->Print(option);
      }
    }

  printf("---------Digits----------------------------------\n");
  printf("   Total digits:       %d\n",(int)fDigits.size());
  if(showdigits)
    for(UInt_t i = 0; i<fDigits.size(); i++) {
      printf("  Digit %3d: %s\n",i,fDigits[i].AsString());
    }
  for(UInt_t j=0; j < fSnarls.size(); j++ ) {
    printf("---------Snarl %d\n",j);
    printf("      Digits:          %d\n",(int)fSnarls[j].size());
    printf("      Trigger time:    %f ns\n",fSnarls[j].GetTriggerTime()/Munits::ns);
    printf("      Trigger Source:  0x%08X\n",fSnarls[j].GetTriggerSource());
    if(fSnarls[j].Size()>0) {
      printf("      Start time:      %f ns\n",fSnarls[j].GetDigit(0)->GetTime()/Munits::ns);
      printf("      Stop  time:      %f ns\n",fSnarls[j].GetDigit(fSnarls[j].Size()-1)->GetTime()/Munits::ns);
    }
    if(showsnarldigits) 
      for(UInt_t i = 0; i<fSnarls[j].size(); i++) {
        printf("  Digit: %s\n",fSnarls[j][i].AsString());
      }
  }

  printf("---------DAQ-------------------------------------\n");
  printf("  TriggerWord: 0x%08X\n",fTotalTriggerWord);
  printf("  TriggerMask: 0x%08X\n",fDaqTriggerMask);
  printf("--------SimEventResult---------------------------\n");
  fResult.Print();
  printf("-------------------------------------------------\n");

  
}


// Simulation Stages:
void SimDetector::DistributePEs( void )
{
  TObjArrayIter itr(fPeList);
  const DigiPE* digipe = NULL;
  double earlyTime = 1e10;
  while( (digipe = dynamic_cast<DigiPE*>(itr.Next()))!=0 ) {
    if(digipe->GetTime()<earlyTime) earlyTime = digipe->GetTime();
  }
  
  // find the timeshift.
  fTimeShift = 0;
  if(earlyTime<0) {
    while((earlyTime+fTimeShift)<0) fTimeShift+=100.*Munits::ns;
    MAXMSG("DetSim",Msg::kWarning,10) << "Negative times in PE list. Shifting all times forward by " << fTimeShift/Munits::ns << " ns" << endl;
  }
  
  //
  // Shift all the times.
  //
  itr.Reset();
  while( (digipe = dynamic_cast<DigiPE*>(itr.Next()))!=0 ) {
    if(fTimeShift>0) {
      // Hack to get times all positive.
      const_cast<DigiPE*>(digipe)->SetTime(digipe->GetTime() + fTimeShift);
    }    

    // Add the hit in.
    this->SortDigiPE( digipe );
  }

}


 
void SimDetector::SortDigiPE( const DigiPE* digipe ) 
{
  if(digipe==0) return;
  PlexPixelSpotId psid;
  PlexStripEndId seid;
  
  const DigiScintHit* hit = digipe->GetHitPointer();
  psid = digipe->GetPixelSpotId();
  if(hit) seid = hit->StripEndId();     // Truth info: convert from DigiScintHit into SEID.
  Double_t time = digipe->GetTime();

  if(!psid.IsValid()) {    
    MSG("DetSim",Msg::kWarning) << "Invalid pixel spot on a DigiPE. Omitting. " << psid.AsString() << "\n";
    MSG("DetSim",Msg::kWarning) << "Time: " << time 
                                << "  DarkNoise: " << ((digipe->IsDarkNoise())?"true":"false")
                                << "  FibreLight: " << ((digipe->IsFibreLight())?"true":"false")
                                << "  ScintHit: \n";
    if(hit) 
      MSG("DetSim",Msg::kWarning) << hit->AsString() << endl;
    
    return;
  }
  if (psid.GetDetector()==Detector::kUnknown) {
    MSG("DetSim",Msg::kWarning) << "Bad pixel spot on a DigiPE (unknown detector). Omitting.\n";
    return;
  } 
  
  // Find the tube.
  SimPmt* pmt = 0;
  int tube = GetTubeIndex(psid);
  SimPmtList::iterator it = fPmts.find(tube);
  pmt = it->second;
  
  // Create the tube if neccessary.
  if(it == fPmts.end()) {
    SimPmt::Pmt_t type = SimPmt::kM16;
    // Determine type of tube (M16 or M64) from the psid.
    if(psid.GetElecType() == ElecType::kQIE) type = SimPmt::kM64;
    
    std::string modelName;
    if( type == SimPmt::kM16 ) modelName = fM16Model;
    else                       modelName = fM64Model;

    pmt = SimPmtMaker::Create(modelName,psid,fContext,fRandom);
    
    if(!pmt) {
      MSG("DetSim",Msg::kError) << "SimDetector: Failed to create new " 
                                << "pmt of type " << modelName << endl;
      return;
    }

    // Call Reset() on the tube to get it up-to-date
    pmt->Reset(fContext);
    
    MSG("DetSim",Msg::kDebug) << "SimDetector: Adding new PMT " 
                              << pmt->GetTubeId().AsString() 
                              << "  type: " << pmt->GetType() 
                              << "  at: " << tube 
                              << endl;
    fPmts[tube] = pmt;
  }
  
  // Add the hit.
  MSG("DetSim",Msg::kVerbose) << "SimDetector::AddDigiPE " 
                              << psid.AsString() << "  "
                              << 1.0 << "pe   " 
                              << time << "ns"
                              << endl;
  
  
  fPmts[tube]->AddDigiPE(digipe);  

  // Add the tube to the 'active' list.
  fActivePmts[tube] = fPmts[tube];
}



void SimDetector::SimulatePmts( void )
{
  SimPmtList::iterator it;
  // Only simulate PMTs that have light. This helps to ensure
  // random-numbers are OK. 
  for(it = fActivePmts.begin(); it!= fActivePmts.end(); it++)  {
    SimPmt* pmt = it->second;
    if(pmt){
      
      MSG("DetSim",Msg::kDebug) << "Simulating PMT "
                                << pmt->GetTubeId().AsString()
                                << endl;
      pmt->SimulatePmt();}
  }
}


void SimDetector::SimulateElectronics( void )
{
  fDigits.clear();

  // Generic FD electronics:
  if(fFarElectronics)
    fFarElectronics->ReadoutDetector(fActivePmts);

  // Now, read them out.
  fDigits.AddList(fFarElectronics->GetDigits());

  // Generic ND electronics:
  if(fNearElectronics)
    fNearElectronics->ReadoutDetector(fActivePmts);

  // Now, read them out.
  fDigits.AddList(fNearElectronics->GetDigits());

  for(UInt_t i = 0; i<fDigits.size(); i++) {
    MSG("DetSim",Msg::kVerbose) 
      << "SimDetector read out digit: " 
      << fDigits[i].AsString() << std::endl;;
  }

}

void SimDetector::SimulateTrigger( void )
{

  
  fTotalTriggerWord = 0;
  fBiggestSnarl = 0;
  fBiggestSnarlSize = -1;

  // Sort the digit list in time for ease in doing the trigger computation.
  fDigits.Sort();

  bool dotrigger = (fApplyTrigger>0);
  static bool donemessage = false;
  if(fApplyTrigger==3)  {
    // Choose an 'intelligent' trigger mode.
    if(fContext.GetDetector()==Detector::kFar) {
      if(!donemessage) MSG("DetSim",Msg::kWarning) 
        << "This is the far detector, so I'm assuming you want DAQ triggers on. Use applytrigger=0 to turn them off." << endl; 
      dotrigger = true;
    }
    else if(fContext.GetDetector()==Detector::kNear) {
       if(!donemessage) MSG("DetSim",Msg::kWarning) 
        << "This is the near detector, so I'm assuming you want DAQ triggers off. Use applytrigger=1 to turn them on." << endl;      
       dotrigger = false;
    }
    else if(fContext.GetDetector()==Detector::kCalDet) {
      if(!donemessage) MSG("DetSim",Msg::kWarning) 
        << "This is the caldet detector, so I'm assuming you want DAQ triggers off. Use applytrigger=1 to turn them on." << endl;      
      dotrigger = false;
    }    
    donemessage = true;
  };

  if(!dotrigger) {
    // Look to see if the all digits pass any triggers, just for 
    // the sake of completeness.
    Int_t triggerWord = 0;
    
    for(int itrig=0;itrig<kSimMaxTriggers;itrig++) {
      int triggerbit = (1<<itrig);
      
      if(fDaqTriggers[itrig]) {
        bool result = fDaqTriggers[itrig]->ApplyTriggerAt(fPlex,
                                                          fDigits,
                                                          0,
                                                          fDigits.size()
                                                          );      
        if(result) triggerWord |= triggerbit;
      } 
    }
    
    // Add trigger bit for SpillTrigger
    if(fContext.GetDetector() == Detector::kFar) {
      triggerWord |= 0x20000; // TRIGGER_BIT_SPILL_IP from mdTriggerCodes.h
    } else
    if(fContext.GetDetector() == Detector::kNear) {
      triggerWord |= 0x10000; // TRIGGER_BIT_SPILL from mdTriggerCodes.h
    }
    
    // Don't bother triggering.. just dump everything into one big snarl.
    double snarltime = 0;
    if(fDigits.size()>0) snarltime = fDigits[0].GetTime();
    SimSnarl snarl(triggerWord,snarltime);
    fTotalTriggerWord = triggerWord;


    snarl.AddList(fDigits);
    
    // Add the snarl.
    fSnarls.push_back(snarl);
    fBiggestSnarlSize = snarl.Size();
    return;
  }

  
  // We now have a sorted list of digits.
  // Go through the list, looking for a burst or a
  // gap, as befits our trigger type.

  UInt_t startdigit = 0;

  while(startdigit<fDigits.size()) {
    int digitsInBurst = 0;
    switch(fDaqTrigTimingMethod) {
    case SimDaqTrigger::kNone:
      digitsInBurst=fDigits.size(); // Always trigger.
      MSG("DetSim",Msg::kDebug) << "Trigger: NULL trigger, triggering all digits." << endl;
      break;
    case SimDaqTrigger::kWindow: // Window trigger.
      digitsInBurst = IsWindowBurst(fDigits,startdigit);
      if(digitsInBurst) MSG("DetSim",Msg::kDebug) << "Trigger: Window trigger fired with " << digitsInBurst << " digits." << endl;
      break;
    case SimDaqTrigger::kGap: // Gap trigger.
      digitsInBurst = IsGapBurst(fDigits,startdigit);
      if(digitsInBurst) MSG("DetSim",Msg::kDebug) << "Trigger: Gap trigger fired with " << digitsInBurst << " digits." << endl;      break;
    default:
      MSG("DetSim",Msg::kError) << "Unknown DAQ trigger window method " << fDaqTrigTimingMethod << endl;
      break;
    }

    if(digitsInBurst>0) { // A burst is starting here.

      // See if this passes any trigger condition:
      Int_t triggerWord = 0;
         
      for(int itrig=0;itrig<kSimMaxTriggers;itrig++) {
        int triggerbit = (1<<itrig);
        
        if(fDaqTriggers[itrig]) {
          bool result = fDaqTriggers[itrig]->ApplyTriggerAt(fPlex,
                                                            fDigits,
                                                            startdigit,
                                                            startdigit+digitsInBurst
                                                            );    
          if(result) triggerWord |= triggerbit;
        }       
      }
      MSG("DetSim",Msg::kDebug) << "Burst digits " << startdigit 
                                << " to " << startdigit+digitsInBurst - 1
                                << " fired triggers: " << triggerWord << endl;

      if(triggerWord != 0){
        // Something triggered.

        fTotalTriggerWord |= triggerWord; // Add to toal.

        // Make a snarl.
        SimSnarl snarl(triggerWord,fDigits[startdigit].GetTime());

        // Collect hits for the trigger.
        Double_t burst_start = fDigits[startdigit].GetTime();
        Double_t burst_stop  = fDigits[startdigit+digitsInBurst-1].GetTime();
        Double_t window_open = burst_start - fDaqTrigReadoutWindowStart;
        Double_t window_close = burst_stop + fDaqTrigReadoutWindowStop;
        
        UInt_t lastdigit = 0;
        for(UInt_t i=0; i<fDigits.size(); i++) {
          if(( fDigits[i].GetTime()>=window_open ) &&
             ( fDigits[i].GetTime()<=window_close ) ) {
            snarl.push_back(fDigits[i]);
            if(i>lastdigit) lastdigit=i;
          }
        }
                
        // Add the snarl.
        
        fSnarls.push_back(snarl);
        if((Int_t)snarl.Size() > fBiggestSnarlSize) {
          fBiggestSnarlSize = snarl.size();;
          fBiggestSnarl = fSnarls.size() -1;
        }

        // Advance the pointer to the end of the burst.
        startdigit += lastdigit;
      }
    }
    startdigit++; // Increment counter by one and look again.
  }
}

void SimDetector::CompileStats() 
{
  fResult.Reset();

  // PMT statistics:
  SimPmtList::iterator it;
  for(it = fActivePmts.begin(); it!= fActivePmts.end(); it++)  {
    SimPmt* pmt = it->second;
    if(pmt){
      fResult.totalPe +=            pmt->GetTotalPe();
      fResult.hitPixels +=          pmt->GetTotalHitPixels(false);
      fResult.hitPixelsWithXtalk += pmt->GetTotalHitPixels(true);
      fResult.totalCharge +=        pmt->GetTotalCharge();
    }
  }

  // Electronics:
  if(fFarElectronics) {
    fResult.digitsAfterFETrigger+=( fFarElectronics->GetDigitsAfterFETrigger() );
    fResult.adcsAfterFETrigger+=(   fFarElectronics-> GetAdcsAfterFETrigger() );
    fResult.digitsAfterSpars+=(     fFarElectronics->GetDigitsAfterSpars() );
    fResult.adcsAfterSpars+=(       fFarElectronics->GetAdcsAfterSpars() );
  }

  // Generic ND electronics:
  if(fNearElectronics) {
    fResult.digitsAfterFETrigger+=( fNearElectronics->GetDigitsAfterFETrigger() );
    fResult.adcsAfterFETrigger+=(   fNearElectronics-> GetAdcsAfterFETrigger() );
    fResult.digitsAfterSpars+=(     fNearElectronics->GetDigitsAfterSpars() );
    fResult.adcsAfterSpars+=(       fNearElectronics->GetAdcsAfterSpars() );
 }

  // Snarls.
  fResult.snarls = GetNumberOfSnarls();
 
  // Create arrays.
  fResult.snarl_digits  = new Int_t[ fResult.snarls ];
  fResult.snarl_trigger = new Int_t[ fResult.snarls ];
  fResult.snarl_adcs    = new Float_t[ fResult.snarls ];
  
  fResult.bigSnarl= GetBiggestSnarl();

  // Loop through snarls. Fill all the stats.
  for(UInt_t isnarl = 0; isnarl < fSnarls.size(); isnarl++) {
    Float_t adctot = 0;
    for(UInt_t i=0; i<fSnarls[isnarl].size(); i++) {
      SimDigit& simdigit = fSnarls[isnarl][i];

      adctot+= simdigit.GetADC() ;
    }
    
    // totals.
    fResult.digitsAfterDaqTrigger+=( fSnarls[isnarl].Size() );
    fResult.adcsAfterDaqTrigger+=( adctot );
    
    fResult.snarl_digits[isnarl] = fSnarls[isnarl].Size();
    fResult.snarl_trigger[isnarl] = fSnarls[isnarl].GetTriggerSource();
    fResult.snarl_adcs[isnarl] = adctot;
  }

  fResult.timeShift = fTimeShift;

}


Int_t SimDetector::IsWindowBurst(SimDigitList& digitList,
                               UInt_t startdigit)
{
  // Assign a winow at startdigit.
  
  // Note: assumes fDaqTrigTimeWindow in TDC units
  // Assumes all electronics identitcal.
  Double_t t_start = digitList[startdigit].GetTime();
  Double_t t_stop  = t_start + fDaqTrigTimeWindow;
  
  Int_t ndigs = 1;
  while(startdigit+ndigs < digitList.GoodSize()) {
    Double_t t = fDigits[startdigit+ndigs-1].GetTime();
    if(t>t_stop) break;
    ndigs++;                    
  }
  return ndigs;
}


Int_t SimDetector::IsGapBurst(SimDigitList& digitList,
                              UInt_t startdigit)
{
  // See if startdigit is the start of a gap trigger.

  // The first digit is always at a gap. Further digits are not at a gap start
  // if they are too far separated from the previous one.
  if(startdigit!=0) {
    Double_t t1 = fDigits[startdigit-1].GetTime();
    Double_t t2 = fDigits[startdigit].GetTime();
    if (t2-t1 < fDaqTrigTimeGap) 
      return 0; // Not the start of a gap.
  }
  
  // Look for gap end.
  Int_t ndigs = 1;
  while(startdigit+ndigs < digitList.GoodSize()) {
    Double_t t1 = fDigits[startdigit+ndigs-1].GetTime();
    Double_t t2 = fDigits[startdigit+ndigs].GetTime();
    if(t2-t1 >= fDaqTrigTimeGap) return ndigs;
    ndigs++;
  }
  return ndigs;
  
}

---