/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes
EsIdealFeeTool Class Reference

#include <EsIdealFeeTool.h>

Inheritance diagram for EsIdealFeeTool:
Inheritance graph
[legend]
Collaboration diagram for EsIdealFeeTool:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 EsIdealFeeTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~EsIdealFeeTool ()
virtual StatusCode generateSignals (DayaBay::ElecPulseCollection *, DayaBay::ElecCrate *)
 This is the extension. Sub classes must provide it.
virtual StatusCode initialize ()
virtual StatusCode finalize ()

Static Public Member Functions

static const InterfaceID & interfaceID ()
 Retrieve interface ID.

Private Member Functions

StatusCode loadResponse ()
StatusCode mapPulsesByChannel (const DayaBay::ElecPulseCollection::PulseContainer &pulses, std::map< DayaBay::ElecChannelId, std::vector< DayaBay::ElecPmtPulse * > > &pulseMap, std::vector< int > &hitCountSlices)
StatusCode updateBoardSignals (DayaBay::ElecFeeCrate *crate)
StatusCode updateEsumSignals (DayaBay::ElecFeeCrate *crate)
virtual StatusCode generateOneChannel (const DayaBay::FeeChannelId &, std::vector< DayaBay::ElecPmtPulse * > &, DayaBay::ElecFeeChannel &, double simTime, const ServiceMode &)
int convertClock (int clock, int inputFrequency, int outputFrequency)
DayaBay::ESumComp::ESumComp_t boardEsumType (int boardId)
void convolve (const DayaBay::AnalogSignal &signal, const DayaBay::AnalogSignal &response, DayaBay::AnalogSignal &output)
void shape (const DayaBay::AnalogSignal &rawSignal, DayaBay::AnalogSignal &shapedSignal, const DayaBay::AnalogSignal &response, std::map< int, int > convolutionWindows)
void digitize (const DayaBay::AnalogSignal &signal, DayaBay::DigitalSignal &output, double range, int bits, double offset)
void discriminate (const DayaBay::AnalogSignal &signal, double threshold, std::vector< int > &output, DayaBay::Threshold::Threshold_t mode)
double pmtPulse (double deltaT, int nPulse)
double overshoot (double deltaT)
double adcOvershoot (double deltaT)
double ringing (double deltaT)
double adcRinging (double deltaT)
double ringingFast (double deltaT)
double ringingSlow (double deltaT, double phase)
double ampSatFactor (int nPulse)
double chargeSatFactor (int nPulse)
double saturationModel (double q, double qSat, double a)
double pulseShaping (double deltaT)
double esumShaping (double deltaT)
double getPulseWidth (int nPulse)
double getPulseForm (int nPulse)
double getPulseShift (int nPulse)
double getOvershootAmp (int nPulse)
double noise ()
void sample (const DayaBay::AnalogSignal &signal, DayaBay::AnalogSignal &output, int inputFrequency, int outputFrequency)

Private Attributes

std::string m_cableSvcName
std::string m_simDataSvcName
ICableSvcm_cableSvc
ISimDataSvcm_simDataSvc
int m_triggerWindowCycles
double m_gainFactor
int m_simFrequency
int m_eSumFrequency
bool m_enableNonlinearity
unsigned int m_linearityThreshold
bool m_enablePretrigger
double m_eSumTriggerThreshold
int m_nHitTriggerThreshold
int m_pretriggerThreshold
bool m_enableNoise
bool m_enableDynamicWaveform
bool m_enableSaturation
bool m_enableRinging
bool m_enableOvershoot
bool m_enableESumTotal
bool m_enableESumH
bool m_enableESumL
bool m_enableFastSimMode
double m_pulseCountSlotWidth
int m_pulseCountWindow
double m_noiseAmp
double m_speAmp
double m_discThreshScale
double m_adcRange
double m_adcBits
double m_adcOffset
DayaBay::AnalogSignal m_pmtPulse
DayaBay::AnalogSignal m_overshoot
DayaBay::AnalogSignal m_ringing
DayaBay::AnalogSignal m_shapedPmtPulse
DayaBay::AnalogSignal m_shapedOvershoot
DayaBay::AnalogSignal m_shapedRinging
DayaBay::AnalogSignal m_esumResponse
DayaBay::AnalogSignal m_shapingResponse
double m_shapingSum
double m_esumShapingSum
double m_pulseMax
Rndm::Numbers m_gauss
DayaBay::AnalogSignal m_rawSignal
DayaBay::AnalogSignal m_esumL
DayaBay::AnalogSignal m_esumH
DayaBay::AnalogSignal m_esumTotal
DayaBay::DigitalSignal m_shapedEsumADC
DayaBay::AnalogSignal m_shapedEsumH
DayaBay::AnalogSignal m_shapedEsumL
DayaBay::AnalogSignal m_shapedEsumTotal
DayaBay::AnalogSignal m_discriminatorSignal
DayaBay::AnalogSignal m_shapedSignal
DayaBay::AnalogSignal m_tdcSignal
DayaBay::AnalogSignal m_adcAnalogSignal
DayaBay::DigitalSignal m_hitSignal
DayaBay::DigitalSignal m_adcHighGain
DayaBay::DigitalSignal m_adcLowGain
DayaBay::AnalogSignal m_energySignal
DayaBay::DigitalSignal m_hitSync
DayaBay::DigitalSignal m_hitHold
std::map< int, int > m_adcConvolutionWindows
std::map< int, int > m_esumConvolutionWindows

Detailed Description

Definition at line 39 of file EsIdealFeeTool.h.


Constructor & Destructor Documentation

EsIdealFeeTool::EsIdealFeeTool ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 11 of file EsIdealFeeTool.cc.

  : GaudiTool(type,name,parent),
    m_cableSvc(0)
{
    declareInterface< IEsFrontEndTool >(this);
    
    declareProperty("CableSvcName",m_cableSvcName="CableSvc",
                    "Name of service to map between detector, hardware, and electronic IDs");
    declareProperty("SimDataSvcName",m_simDataSvcName="SimDataSvc",
                    "Name of service to provide FEE channel properties for simulation");
    declareProperty("SimFrequency", m_simFrequency = DayaBay::TdcFrequencyHz,
                    "Simulation Frequency");
    declareProperty("ESumFrequency", m_eSumFrequency = DayaBay::EsumFrequencyHz,
                    "Esum Frequency");
    declareProperty("TriggerWindowCycles", 
                    m_triggerWindowCycles=DayaBay::TriggerWindowCycles,
                    "Number of trigger window ccyles");
    //  Nov. 30 Changed the gain to 1e7. -- Yasu.
    declareProperty("GainFactor", 
                    m_gainFactor=1.0,
                    "Mean PMT gain (in units of 10^7)");
    declareProperty("NoiseAmp",m_noiseAmp=0.5e-3*Gaudi::Units::volt,
                    "Voltage amplitude of noise");
    declareProperty("SpeAmp",m_speAmp=0.0035,
                    "Spe pulse height in V at 1e7 gain");
    //  Nov. 30: Changed to the discriminater threshold to 0.15 p.e., which corresponds to FEE threshold of 300, accoring to doc:6963
    // The original threshould is 0.25 p.e. defined in $DATASVCROOT/share/feeDataTable.txt 
    // So apply factor of 0.6 (= 0.15/0.25)    -- Yasu
    // change it from 0.6 to 0.5, corresponds to FEE threshold of 250, (2012.02.16) by Guofu
    declareProperty("DiscThreshScale", m_discThreshScale=0.5,
                    "Factor to scale all channels' discriminator threshold");
    declareProperty("EnableNonlinearity", m_enableNonlinearity=true, 
                    "Turns on/off all nonlinear effects");
    declareProperty("EnableNoise",m_enableNoise=false,"Turn on/off noise individually");
    declareProperty("EnableDynamicWaveform", m_enableDynamicWaveform=true, 
                    "Turns on/off dynamic waveform individually");
    declareProperty("EnableRinging", m_enableRinging=true, 
                    "Turns on/off ringing individually");
    declareProperty("EnableOvershoot", m_enableOvershoot=true, 
                    "Turns on/off overshoot individually");
    declareProperty("EnableSaturation", m_enableSaturation=true, 
                    "Turns on/off saturation individually");
    declareProperty("EnableESumTotal", m_enableESumTotal=true, 
                    "Turns on/off Total ESum simulation");
    declareProperty("EnableESumH", m_enableESumH=true, "Turns on/off upper ESum simulation");
    declareProperty("EnableESumL", m_enableESumL=true, "Turns on/off lower ESum simulation");
    declareProperty("EnableFastSimMode", m_enableFastSimMode=false, 
                    "Turns on/off fast simulation mode to save CPU time");
    declareProperty("PulseCountSlotWidth", m_pulseCountSlotWidth=10e-9*Gaudi::Units::second,
                    "pulse count time slot size");
    declareProperty("PulseCountWindow", m_pulseCountWindow=1, 
                    "Number of counted slots before and after a pulse");
    declareProperty("LinearityThreshold", m_linearityThreshold=20, 
                    "Threshold in P.E. for activating non-linear pulse shape");
    declareProperty("EnablePretrigger", m_enablePretrigger=true, 
                    "Turns on/off pretrigger");
    declareProperty("TotalESumTriggerThreshold", 
                    m_eSumTriggerThreshold=DayaBay::Trigger::kADESumThreshold, 
                    "Total ESum trigger threshold as set in TrigSim");
    declareProperty("NHitTriggerThreshold", 
                    m_nHitTriggerThreshold=DayaBay::Trigger::kADthreshold, 
                    "nHit trigger threshold as set in TrigSim");
    declareProperty("ADCOffset", m_adcOffset = 500, "Esum Frequency");
    declareProperty("ADCRange", m_adcRange = 1, "Esum Range");
    declareProperty("ADCBits", m_adcBits = 10, "Esum Bits");
}
EsIdealFeeTool::~EsIdealFeeTool ( ) [virtual]

Definition at line 80 of file EsIdealFeeTool.cc.

{}

Member Function Documentation

StatusCode EsIdealFeeTool::generateSignals ( DayaBay::ElecPulseCollection ,
DayaBay::ElecCrate  
) [virtual]

This is the extension. Sub classes must provide it.

Implements IEsFrontEndTool.

Definition at line 157 of file EsIdealFeeTool.cc.

{    
  // Ensure that crate is a FeeCrate
  DayaBay::ElecFeeCrate* crate = 0;
  crate = dynamic_cast<DayaBay::ElecFeeCrate*>(elecCrate);
  if(crate == NULL){
    error() << "generateSignals(): Crate is not a FEE Crate." << endreq;
    return StatusCode::FAILURE;
  }
  
  // Context, ServiceMode for this data
  Context context(pulses->detector().site(), SimFlag::kMC, 
          pulses->header()->header()->timeStamp(),
          pulses->detector().detectorId());
  int task = 0;
  ServiceMode svcMode(context, task);

  // Initialize crate if necessary
  if(crate->channelData().size() == 0){
    // First, get list of all connected FEE channels
    const std::vector<DayaBay::DetectorSensor>& pmtList 
      = m_cableSvc->sensors( svcMode );
    int nPmts = pmtList.size();
    for(int pmtIdx = 0; pmtIdx < nPmts; pmtIdx++){
      // Add channel to crate
      DayaBay::FeeChannelId channelId(m_cableSvc->elecChannelId(pmtList[pmtIdx], svcMode).fullPackedData());
      crate->addChannel(channelId);
      verbose() << "Adding channel to crate: " << channelId << endreq;
    }
  }

  // Prepare time range of the simulation
  TimeStamp earliest = pulses->header()->header()->earliest();
  TimeStamp latest = pulses->header()->header()->latest();
  // The time window for the crate signals
  double simTime = latest - earliest;
  verbose() << "Simulation time is " << simTime << endreq;
  // The number of samples in time given the simulation frequency
  int simSamples = int(simTime * m_simFrequency);
  verbose() << "Simulation has " << simSamples << " samples at frequency "
     << m_simFrequency << " Hz" << endreq;
  
  // Compute time windows that will possibly be readout
  int nTimeSlices = int(simTime / 150e-9 + 0.5);
  std::vector<int> hitCountSlices(nTimeSlices,0);
  for ( int i = 0; i < nTimeSlices; ++i ) hitCountSlices[i] = 0;

  // Organize Pulses by Channel
  std::map<DayaBay::ElecChannelId, std::vector<DayaBay::ElecPmtPulse*> > pulseMap;
  StatusCode sc;
  sc = mapPulsesByChannel( pulses->pulses(), pulseMap, hitCountSlices );      
  if(sc.isFailure()) return sc;
  debug() << "Mapped Pulses" << endreq;

  // Pretrigger:
  // Compute map of time windows that will possibly read out  
  m_adcConvolutionWindows.clear();
  m_esumConvolutionWindows.clear();
  int adcConvStart, adcConvEnd;
  int esumConvStart, esumConvEnd;
  std::map<int,int>::reverse_iterator sliceIt;
  for(int i = 1; i < nTimeSlices; i++){
    if( hitCountSlices[i] + hitCountSlices[i-1] < m_pretriggerThreshold ) continue;
    if (i - 5 < 0) adcConvStart = 0;
    else adcConvStart = (i - 5);
    if (i - 1 < 0) esumConvStart = 0;
    else esumConvStart = (i - 1);
    
    if (i + 5 > nTimeSlices) adcConvEnd = nTimeSlices;
    else adcConvEnd = (i + 5);
    if (i + 1 > nTimeSlices) esumConvEnd = nTimeSlices;
    else esumConvEnd = (i + 1);
    
    if( m_adcConvolutionWindows.empty() ) {
      m_adcConvolutionWindows[adcConvStart] = adcConvEnd;
      m_esumConvolutionWindows[esumConvStart] = esumConvEnd;
      continue;
    }
    sliceIt = m_adcConvolutionWindows.rbegin();
    if( (adcConvStart - sliceIt->second > 2) ) 
      m_adcConvolutionWindows[adcConvStart ] = adcConvEnd;
    else
      m_adcConvolutionWindows[sliceIt->first] = adcConvEnd;
      
    sliceIt = m_esumConvolutionWindows.rbegin();
    if( (esumConvStart - sliceIt->second > 2) ) 
      m_esumConvolutionWindows[esumConvStart ] = esumConvEnd;
    else                                         
      m_esumConvolutionWindows[sliceIt->first] = esumConvEnd;
  }
  debug() << "Number of time windows to be possibly read out: " << m_adcConvolutionWindows.size() << endreq;

  // If no time windows is considered of possible interest, skip the entire event
  if(m_adcConvolutionWindows.size()==0 && m_enablePretrigger) {
    debug() << "Event below pretrigger threshold, no waveform simulation" << endreq;
    return StatusCode::SUCCESS;
  }
  
  // Loop over channels, and fill with the simulated signals
  const DayaBay::ElecFeeCrate::ChannelData& channelData = crate->channelData();
  DayaBay::ElecFeeCrate::ChannelData::const_iterator channelIter, 
    done = channelData.end();  
  for(channelIter = channelData.begin(); channelIter != done; ++channelIter){
    DayaBay::FeeChannelId channelId = channelIter->first;
    DayaBay::ElecFeeChannel& channel = crate->channel(channelId);

    //Fill each channel with signals
    sc = generateOneChannel(channelId, pulseMap[channelId], channel, 
                simTime, svcMode);
    if( !sc.isSuccess() ) return sc;

  } // End loop over channels

  StatusCode status = updateBoardSignals(crate);
  if (status.isFailure()) return status;
  status = updateEsumSignals(crate);
  if (status.isFailure()) return status;  

  verbose() << "Crate\n" << *crate << endreq;

  return StatusCode::SUCCESS;
}
StatusCode EsIdealFeeTool::initialize ( ) [virtual]

Definition at line 82 of file EsIdealFeeTool.cc.

{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);

  // Get PMT simulation input data service
  m_simDataSvc = svc<ISimDataSvc>(m_simDataSvcName,true);

  StatusCode sc = loadResponse();
  if(sc.isFailure()) return sc;

  // Random number service
  IRndmGenSvc *rgs = 0;
  if (service("RndmGenSvc",rgs,true).isFailure()) {
    fatal() << "Failed to get random service" << endreq;
    return StatusCode::FAILURE;        
  }
  
  sc = m_gauss.initialize(rgs, Rndm::Gauss(0,1));
  if (sc.isFailure()) {
    fatal() << "Failed to initialize gaussian random numbers" << endreq;
    return StatusCode::FAILURE;
  }
  // Set up pretrigger threshold
  int preThr = int(0.8 * m_nHitTriggerThreshold + 0.5);
  int esumPreThr = int (0.8 * 1500 * m_eSumTriggerThreshold/Gaudi::Units::volt + 0.5);
  if(preThr > esumPreThr) preThr = esumPreThr;
  m_pretriggerThreshold = preThr;
  if( m_enablePretrigger ){
    info() << "Pretrigger threshold is " << m_pretriggerThreshold << " pe " << endreq;
  }
  else {
    info() << "Pretrigger disabled" << endreq;
    info() << "Convolution threshold is " << m_pretriggerThreshold << " pe " << endreq;
  }
  
  if( m_enableFastSimMode){
    info() << "Fast electronics simulation mode enabled: Note that Upper and Lower ESum trigger do not yield meaningful results" << endreq;
    m_enableDynamicWaveform = false;
    m_enableESumL = false;
    m_enableESumL = false;
  }
        
  if( !m_enableOvershoot ) 
    info() << "Overshoot disabled" << endreq; 
  else 
    debug() << "Overshoot enabled" << endreq;
    
  if( !m_enableRinging ) 
    info() << "Ringing disabled" << endreq; 
  else 
    debug() << "Ringing enabled" << endreq;
    
  if( !m_enableSaturation ) 
    info() << "Saturation disabled" << endreq; 
  else 
    debug() << "Saturation enabled" << endreq;
    
  if( !m_enableDynamicWaveform ) {
    info() << "Dynamic waveform disabled" << endreq; 
    m_linearityThreshold = 1e6;
  }
  else 
    info() << "Dynamic waveform enabled for channels with more than " 
         << m_linearityThreshold << " hits " << endreq;

  
  return StatusCode::SUCCESS;
}
StatusCode EsIdealFeeTool::finalize ( ) [virtual]

Definition at line 152 of file EsIdealFeeTool.cc.

{
  return StatusCode::SUCCESS;
}
StatusCode EsIdealFeeTool::loadResponse ( ) [private]

Definition at line 503 of file EsIdealFeeTool.cc.

                                       { 
  double dT_seconds = (1. / m_simFrequency);
  double dT_units = dT_seconds*CLHEP::second;
  int nShapingSamples = int(1.2*DayaBay::preTimeTolerance/dT_units);
  int nPulseSamples = int(1.2*DayaBay::preTimeTolerance/dT_units);
  int nOvershootSamples = int(16000/dT_units);
  DayaBay::AnalogSignal overshootAdc;
  DayaBay::AnalogSignal ringingAdc;
  m_pmtPulse.resize(nPulseSamples);
  m_overshoot.resize(nOvershootSamples);
  m_ringing.resize(nOvershootSamples);
  m_shapingResponse.resize(nShapingSamples);
  m_esumResponse.resize(nPulseSamples);
  overshootAdc.resize(nOvershootSamples);
  ringingAdc.resize(nOvershootSamples);
  m_shapingSum = 0;
  m_esumShapingSum = 0;
  m_pulseMax = 0;

  for (int i=0; i<nShapingSamples; i++) {
    m_shapingResponse[i] = pulseShaping(i*dT_seconds);
    m_shapingSum += m_shapingResponse[i];
  }
  
  for (int i=0; i<nPulseSamples; i++){
    m_esumResponse[i] = esumShaping(i*dT_units/CLHEP::nanosecond);
    m_esumShapingSum += m_esumResponse[i];
  }

  for (int i=0; i<nPulseSamples; i++) {
    m_pmtPulse[i] = pmtPulse(i*dT_seconds,1);
  }
  
  for (int i=0; i<nOvershootSamples; i++) {
    m_overshoot[i] = 0;
    m_ringing[i] = 0;
    // Store overshoot/ringing for raw signal
    if(m_enableOvershoot) m_overshoot[i] += overshoot(i*dT_seconds);
    if(m_enableRinging)   m_ringing[i]   += ringing(i*dT_seconds);
    // Overshoot/ringing in ADC signal after FEE shaping has a different shape
    overshootAdc[i] = 0;
    ringingAdc[i]   = 0;
    if(m_enableOvershoot) overshootAdc[i]+= adcOvershoot(i*dT_seconds);
    if(m_enableRinging)   ringingAdc[i]  += adcRinging(i*dT_seconds);
  }

  // Convolve the pmt pulse, overshoot, and ringing 
  // with the CR-(RC)^4 shaping response
  convolve( m_pmtPulse, m_shapingResponse, m_shapedPmtPulse );
  m_shapedPmtPulse.resize( m_pmtPulse.size() );
  convolve( overshootAdc, m_shapingResponse, m_shapedOvershoot );
  m_shapedOvershoot.resize( m_overshoot.size() );
  convolve( ringingAdc, m_shapingResponse, m_shapedRinging );
  m_shapedRinging.resize( m_ringing.size() );
  // Scale the resulting convolution to preserve the pulse area
  for(unsigned int i=0; i < m_pmtPulse.size(); i++) 
    m_shapedPmtPulse[i] /= m_shapingSum;
  for(unsigned int i=0; i < m_overshoot.size(); i++) 
  {
    m_shapedOvershoot[i] /= m_shapingSum;
    m_shapedRinging[i]  /= m_shapingSum;
  }
  // Dump ideal pulse shapes
  {
    std::stringstream output;
    output << "pmtPulse = [ ";
    for(int i=0; i<nPulseSamples; i++){
      if(i!=0) output << ", ";
      output << m_pmtPulse[i];
    } 
    output << " ]";
    verbose() << output.str() << endreq;
  }
  {
    std::stringstream output;
    output << "overshoot = [ ";
    for(int i=0; i<nOvershootSamples; i++){
      if(i!=0) output << ", ";
      output << m_overshoot[i];
    } 
    output << " ]";
    verbose() << output.str() << endreq;
  }
  
  {
    std::stringstream output;
    output << "Esum Shaping Pulse = [ ";
    for(int i=0; i<nPulseSamples; i++){
        if(i!=0) output << ", ";
    output << m_esumResponse[i];
    }
    output << " ]";
    verbose() << output.str() << endreq;
  }
  
  {
    std::stringstream output;
    output << "shapedPmtPulse = [ ";
    for(int i=0; i<nPulseSamples; i++){
      if(i!=0) output << ", ";
      output << m_shapedPmtPulse[i];
    } 
    output << " ]";
    verbose() << output.str() << endreq;
  }
  {
    std::stringstream output;
    output << "shapedOvershoot = [ ";
    for(int i=0; i<nOvershootSamples; i++){
      if(i!=0) output << ", ";
      output << m_shapedOvershoot[i];
    } 
    output << " ]";
    verbose() << output.str() << endreq;
  }
  return StatusCode::SUCCESS;
}
StatusCode EsIdealFeeTool::mapPulsesByChannel ( const DayaBay::ElecPulseCollection::PulseContainer pulses,
std::map< DayaBay::ElecChannelId, std::vector< DayaBay::ElecPmtPulse * > > &  pulseMap,
std::vector< int > &  hitCountSlices 
) [private]

Definition at line 621 of file EsIdealFeeTool.cc.

{
  debug() << "Processing " << pulses.size() << " pmt pulses." << endreq;
  DayaBay::ElecPulseCollection::PulseContainer::const_iterator it, 
    pulseDone = pulses.end();
  int timeSliceNo;
  // loop over pulses in collection
  for (it=pulses.begin(); it != pulseDone; ++it) {
    DayaBay::ElecPmtPulse* pulse = dynamic_cast<DayaBay::ElecPmtPulse*>(*it);
    if(pulse == NULL){
      // Catch invalid pulses
      error() << "mapPulsesByChannel(): bad pulse." << endreq;
      return StatusCode::FAILURE;
    }
    (pulseMap[pulse->channelId()]).push_back(pulse);
    timeSliceNo = int( pulse->time()/150. - 0.5 );
    hitCountSlices[timeSliceNo]+=int(pulse->amplitude()+0.5);
  }
  return StatusCode::SUCCESS;
}
StatusCode EsIdealFeeTool::updateBoardSignals ( DayaBay::ElecFeeCrate crate) [private]

Definition at line 645 of file EsIdealFeeTool.cc.

{ 
  // Update the board Nhit and Esum info, according to the current
  // channel info
  
  DayaBay::ElecFeeCrate::DigitalMap& nhitMap = crate->nhit();
  DayaBay::ElecFeeCrate::AnalogMap& esumMap = crate->esum();

  // Loop over channels, and add data to the board-level signals
  const DayaBay::ElecFeeCrate::ChannelData& channelData = crate->channelData();
  DayaBay::ElecFeeCrate::ChannelData::const_iterator channelIter, 
    done = channelData.end();
  for(channelIter = channelData.begin(); channelIter != done; ++channelIter){
    DayaBay::FeeChannelId channelId = channelIter->first;
    const DayaBay::ElecFeeChannel& channel = channelIter->second;
    const DayaBay::DigitalSignal& hit = channel.hit();
    const DayaBay::AnalogSignal& energy = channel.energy();
    DayaBay::DigitalSignal& boardNhit = nhitMap[channelId.boardId()];
    DayaBay::AnalogSignal& boardEsum = esumMap[channelId.boardId()];
    
    unsigned int hitSamples = hit.size();

    // Create hitSync to mimic the syncronized discriminator output
    //DayaBay::DigitalSignal hitSync( hitSamples,0);
    if( m_hitSync.size() != hitSamples ) m_hitSync.resize( hitSamples );
    const int* hitStart = &hit[0];
    int* hitSyncStart = &m_hitSync[0];
    *hitSyncStart = 0;  // Set first sample to zero
    for (unsigned int i=1; i< hitSamples; i++){
      int* curHitSync = hitSyncStart + i;
      int* lastHitSync = curHitSync - 1;
      const int* curHit = hitStart + i;
      const int* lastHit = curHit - 1;
      if( *lastHitSync==0 && *lastHit==0 && *curHit==1 ) *curHitSync = 1;
      else *curHitSync = 0;
    }
    
    // The logic pulse from the channel level hit is extended for
    // multiple clock cycles (triggerWindowCycles) before going to the
    // trigger.
    //DayaBay::DigitalSignal hitHold( hitSamples , 0 );
    if( m_hitHold.size() != hitSamples ) m_hitHold.resize( hitSamples );
    int* hitHoldStart = &m_hitHold[0];
    for (unsigned int i=0; i< hitSamples; i++) *(hitHoldStart + i) = 0;
    
    if( boardNhit.size() != hitSamples ){
      boardNhit.resize( hitSamples );
      int* boardNhitStart = &boardNhit[0];
      for (unsigned int i=0; i< hitSamples; i++) *(boardNhitStart + i) = 0;
    }
      
    for( unsigned int idx = 0; idx < hitSamples; idx++ ){
      unsigned int tempIdx = 0;
      if( *(hitSyncStart+idx) != 0 ){
    for( int holdIdx = 0; holdIdx < m_triggerWindowCycles; holdIdx++ ){
      unsigned int currentIdx = idx + holdIdx;
      if( currentIdx < hitSamples ) *(hitHoldStart + currentIdx) = 1;
      tempIdx++;
    }
    if(tempIdx !=0) idx = idx +tempIdx-1 ;
      } // End hit
    }

    {
      int* boardNhitStart = &boardNhit[0];
      for(unsigned int idx = 0; idx < hitSamples; idx++)
    *(boardNhitStart + idx) += *(hitHoldStart + idx);
    }
    verbose() << "boardNhit " << boardNhit << endreq;
    
    // Add channel energy to board esum
    unsigned int esumSamples = energy.size();
    if( boardEsum.size() != esumSamples ){
      boardEsum.resize( esumSamples );
      double* boardEsumStart = &boardEsum[0];
      for (unsigned int i=0; i< esumSamples; i++) *(boardEsumStart + i) = 0;
    }
    const double* energyStart = &energy[0];
    double* boardEsumStart = &boardEsum[0];
    for( unsigned int idx = 0; idx < esumSamples; idx++ ){
      *(boardEsumStart + idx) += *(energyStart + idx);
    }
  } // End loop over channels

  return StatusCode::SUCCESS;
}
StatusCode EsIdealFeeTool::updateEsumSignals ( DayaBay::ElecFeeCrate crate) [private]

Definition at line 732 of file EsIdealFeeTool.cc.

                                                                      {
  // Loop over all boards in map
  DayaBay::ElecFeeCrate::AnalogMap& esumMap = crate->esum();
  DayaBay::ElecFeeCrate::AnalogMap::iterator boardIt = esumMap.begin();
  
  // Clear transient stores from previous Execution cycles
  m_esumH.clear();
  m_esumL.clear();
  m_esumTotal.clear();
  
  for (;boardIt != esumMap.end(); ++boardIt)
  {
      DayaBay::FeeChannelId boardId = boardIt->first;
      DayaBay::AnalogSignal& boardEnergy = boardIt->second;
      int boardNum = boardId.board();

      if( m_esumL.size() != boardEnergy.size() ){
          m_esumL.resize( boardEnergy.size() );
      }
      
      if( m_esumH.size() != boardEnergy.size() ){
        m_esumH.resize( boardEnergy.size() );
      }
      
      if( m_esumTotal.size() != boardEnergy.size() ){
        m_esumTotal.resize( boardEnergy.size() );
      }

      debug() << "board number " << boardNum << endreq;
      debug() << "Low Esum size " << m_esumL.size() << endreq;
      debug() << "board Energy size " << boardEnergy.size() << endreq;
      
      int boardIdx = boardEnergy.size();
      double* boardEsumStart = &boardEnergy[0];
      double* esumHStart = &m_esumH[0];
      double* esumLStart = &m_esumL[0];
      double* esumTtart  = &m_esumTotal[0];
      
      for( int i = 0; i < boardIdx; i++)
      {
          // update Upper Sum and total sum
          if (boardEsumType(boardNum) == DayaBay::ESumComp::kESumHigh){
              *(esumHStart+i) += *(boardEsumStart+i);
              *(esumTtart+i) += *(boardEsumStart+i);
              //m_esumH[i] += boardEnergy[i];
              //m_esumTotal[i] += boardEnergy[i];
          }
          
          // update Lower Sum and total sum
          if (boardEsumType(boardNum) == DayaBay::ESumComp::kESumLow){
              *(esumLStart+i) += *(boardEsumStart+i);
              *(esumTtart+i) += *(boardEsumStart+i);
              //m_esumL[i] += boardEnergy[i];
              //m_esumTotal[i] += boardEnergy[i];
          }
      }

  
  }// END LOOP over boards
    
  
  verbose() << "Lower ESum " << m_esumL << endreq;
  verbose() << "Total ESum " << m_esumTotal << endreq;
    
  // Shape Upper Sum
  if(m_enableESumH) {
    verbose() << "Upper ESum " << m_esumH << endreq;
    shape( m_esumH, m_shapedEsumH, m_esumResponse, m_esumConvolutionWindows );
    m_shapedEsumH.resize( m_esumH.size() );
    // NOTE: They should all be the same lenght so if we need to speed things up
    //        we could resacle all 3 in the same loop...
    
    // rescale after convolution
    double* esumHStart = &m_shapedEsumH[0];
    
    for(unsigned int i=0; i < m_shapedEsumH.size(); i++) 
    {  
      *(esumHStart+i) /= m_esumShapingSum;
      //m_shapedEsumH[i] /= m_esumShapingSum;  
    }
  
    // sample down to appropriate frequency and add to header 
    sample( m_shapedEsumH, crate->esumUpper(),
             m_simFrequency, m_eSumFrequency );
    
    verbose() << "Shaped upper sum " << crate->esumUpper() << endreq;
  }
  else {
    m_shapedEsumH.resize( m_esumH.size() );
    crate->setEsumUpper(m_shapedEsumH);
  }
  // Shape Lower Sum
  if(m_enableESumL) {
    shape( m_esumL, m_shapedEsumL, m_esumResponse, m_esumConvolutionWindows );
    m_shapedEsumL.resize( m_esumL.size() );
  
    // rescale after convolution
    double* esumLStart = &m_shapedEsumL[0];
    
    for(unsigned int i=0; i < m_shapedEsumL.size(); i++)
    {
        *(esumLStart+i) /= m_esumShapingSum;
        //m_shapedEsumL[i] /= m_esumShapingSum;
    }
    
    // sample down to appropriate frequency and add to header  
    sample( m_shapedEsumL, crate->esumLower(),
             m_simFrequency, m_eSumFrequency );
      
    verbose() << "Shaped lower sum " << crate->esumLower() << endreq;
  }
  else {
    m_shapedEsumL.resize( m_esumL.size() );
    crate->setEsumLower(m_shapedEsumL);
  }
  // Shape Total Sum 
  if(m_enableESumTotal) {
    shape( m_esumTotal, m_shapedEsumTotal, m_esumResponse, m_esumConvolutionWindows );
    m_shapedEsumTotal.resize( m_esumTotal.size() );
    
    // rescale after convolution
    double* esumTstart  = &m_shapedEsumTotal[0];
     
    for(unsigned int i=0; i < m_shapedEsumTotal.size(); i++) 
    {
       *(esumTstart+i) /= m_esumShapingSum;
       //m_shapedEsumTotal[i] /= m_esumShapingSum; 
    }  
  
    // sample down to appropriate frequency and add to header  
    sample( m_shapedEsumTotal, crate->esumTotal(),
             m_simFrequency, m_eSumFrequency );
        
    verbose() << "Shaped total sum " << crate->esumTotal() << endreq;
  }
  else {
    m_shapedEsumTotal.resize( m_esumTotal.size() );
    crate->setEsumTotal(m_shapedEsumTotal);
  }
  //Digitize the shaped total sum
  //FIXME: This needs to be fixed (edraeger@iit.edu)
  //         set as properties and/or get from service
  digitize(crate->esumTotal(), crate->esumADC(), m_adcRange, m_adcBits, m_adcOffset);
  
  verbose() << "Digitized signal " << crate->esumADC() << endreq;
  
  return StatusCode::SUCCESS;
}
StatusCode EsIdealFeeTool::generateOneChannel ( const DayaBay::FeeChannelId channelId,
std::vector< DayaBay::ElecPmtPulse * > &  channelPulses,
DayaBay::ElecFeeChannel channel,
double  simTime,
const ServiceMode svcMode 
) [private, virtual]

Definition at line 281 of file EsIdealFeeTool.cc.

{
  // Generate the FEE signals caused by the PMT pulses for one channel
  // Get the simulation properties for this channel
  const DayaBay::FeeSimData* feeSimData = 
    m_simDataSvc->feeSimData(channelId, svcMode);
  if(!feeSimData){
    error() << "No Simulation input properties for FEE channel: " 
        << channelId << endreq;
    return StatusCode::FAILURE;
  }

  // The number of samples in time given the simulation frequency
  unsigned int simSamples = int(simTime * m_simFrequency);

  // Prepare Raw Signal
  if(m_rawSignal.size() != simSamples) m_rawSignal.resize( simSamples );
  if(m_shapedSignal.size() != simSamples) m_shapedSignal.resize( simSamples );
  double* rawStart = &m_rawSignal[0];
  double* shapedStart = &m_shapedSignal[0];
  for( unsigned int sigIdx = 0; sigIdx!=simSamples; sigIdx++){
    *(rawStart + sigIdx) = 0;
    *(shapedStart + sigIdx) = 0;
  } 
  
  // Add noise to raw signal if requested
  if( m_enableNoise ){
    for(unsigned int index=0; index < simSamples; index++){
      *(rawStart + index) += noise();
    }
  }
  int channelPulsesN = channelPulses.size();
  // Pulses on this channel?
  if( channelPulsesN > 0 ){
    verbose() << "Channel " << channelId << " has " << channelPulses.size() 
          << " pulses." << endreq;

    // Prepare time slots for pulse counting
    int numPulseTimeSlots = int(simTime*Gaudi::Units::second
                /m_pulseCountSlotWidth) + 1;
    verbose() << "Number of time slots for pulse counting = " 
      << numPulseTimeSlots << endreq; 
    std::vector<int> pulseTimeSlots(numPulseTimeSlots);

    // Fill pulse count time slots
    // Count the number of pulses for each time slot to model nonlinearity
    for (int i = 0; i < numPulseTimeSlots; i++ ) pulseTimeSlots[i] = 0;
    std::vector<DayaBay::ElecPmtPulse*>::const_iterator pulseIter, pulseDone = channelPulses.end();
    for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
      DayaBay::ElecPmtPulse* pulse = *pulseIter;
      int timeSlot = int(pulse->time() * Gaudi::Units::nanosecond 
                 / m_pulseCountSlotWidth);
      int pulseNo = int(pulse->amplitude()+0.5);
      if(timeSlot>=numPulseTimeSlots) continue;
      pulseTimeSlots[timeSlot]+=pulseNo;
      verbose() << "Added " << pulseNo << " pulses in time slot " << timeSlot << endreq;
    }

    double* pmtPulseStart = &m_pmtPulse[0];
    double* overshootStart = &m_overshoot[0];
    double* ringingStart = &m_ringing[0];
    double* shapedPmtPulseStart = &m_shapedPmtPulse[0];
    double* shapedOvershootStart = &m_shapedOvershoot[0];
    double* shapedRingingStart = &m_shapedRinging[0];
    
    unsigned int nPulseSamples = m_pmtPulse.size();
    unsigned int nOvershootSamples = m_overshoot.size();
    
    // Loop over pulses, add each to raw signal
    for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
      DayaBay::ElecPmtPulse* pulse = *pulseIter;
      int tOffset = int(pulse->time()*1e-9 * m_simFrequency);
      float amplitude = pulse->amplitude() * m_speAmp * m_gainFactor;
      unsigned int nPulse=0;
      unsigned int nCoincPulse=0;
      float satAmp = amplitude;
      float satCharge = amplitude;
      // Count the total number of pulses within a nearby time window
      // This number is used to determine the nonlinearity
      if(channelPulsesN>5){
        int timeSlot = int(pulse->time() * Gaudi::Units::nanosecond 
                 / m_pulseCountSlotWidth);
        nCoincPulse = pulseTimeSlots[timeSlot];
        for (int iSlot = timeSlot - m_pulseCountWindow; 
             iSlot <= timeSlot + m_pulseCountWindow; iSlot++){
          if(iSlot>=0 && iSlot<numPulseTimeSlots) 
            nPulse += pulseTimeSlots[iSlot];
        }
        // Get saturated amplitudes
        satAmp *= ampSatFactor(nPulse);
        satCharge *= chargeSatFactor(nPulse);
      }else{
        nCoincPulse =1;
        nPulse = 1;
      }
      
      verbose() << "Number of pulses in time window: " << nPulse << endreq;
      // Cut off overshoot/ringing after 1.5 us for small pe numbers
      unsigned int nSamples = int(3.5 * nPulseSamples);
      // Skip overshoot and ringing if there is only 1 pe 
      // For large pe numbers, add separately per cluster
      if(channelPulsesN ==1 || nCoincPulse > 5
         || (!m_enableOvershoot && !m_enableRinging))
        nSamples = nPulseSamples;
            
      // Now add pulses to the raw and shaped signal
      for(unsigned int index = 0; index < nSamples; index++){
        unsigned int sampleIdx = tOffset + index;
        double idxTime = index * (1. / m_simFrequency);
        if(sampleIdx>0 && sampleIdx<simSamples){
          if(index<nPulseSamples) {
            if(nPulse>m_linearityThreshold)
              *(rawStart + sampleIdx) -= satAmp * pmtPulse(idxTime, nPulse);
            else
              *(rawStart + sampleIdx) -= satCharge * (*(pmtPulseStart + index));
            *(shapedStart + sampleIdx) -= satCharge * (*(shapedPmtPulseStart+index));
          }
          // Add overshoot/ringing
          if((m_enableOvershoot || m_enableRinging) && nCoincPulse < 6){
            *(rawStart + sampleIdx) -= satCharge * (*(overshootStart+index));
            *(rawStart + sampleIdx) -= satAmp * (*(ringingStart+index));
            *(shapedStart + sampleIdx) -= satCharge * (*(shapedOvershootStart+index));
            *(shapedStart + sampleIdx) -= satAmp * (*(shapedRingingStart+index));
          }
        }
      }
    } // End loop over pulses
    
    // Add overshoot/ringing for pe cluster
    for (int timeSlot = 0; timeSlot < numPulseTimeSlots; timeSlot++ ) {
      if (pulseTimeSlots[timeSlot] > 5) {
        
        int nPulse=0;
        for (int iSlot = timeSlot - m_pulseCountWindow; 
             iSlot <= timeSlot + m_pulseCountWindow; iSlot++){
          if(iSlot>=0 && iSlot<numPulseTimeSlots) nPulse += pulseTimeSlots[iSlot];
        }
        float satCharge = pulseTimeSlots[timeSlot] * chargeSatFactor(nPulse) * m_speAmp * m_gainFactor;
        float satAmp = pulseTimeSlots[timeSlot] * ampSatFactor(nPulse) * m_speAmp * m_gainFactor;
        int osStartIdx = int((timeSlot) * m_pulseCountSlotWidth * m_simFrequency * 1e-9);
        unsigned int osEnd = nOvershootSamples;
        if( osEnd + osStartIdx > simSamples ) 
          osEnd = simSamples - osStartIdx;
        for (unsigned int index=0; index<osEnd; index++) {
          int sampleIdx = osStartIdx + index;
          *(rawStart + sampleIdx)    -= satCharge * (*(overshootStart+index));
          *(rawStart + sampleIdx)    -= satAmp    * (*(ringingStart+index));
          *(shapedStart + sampleIdx) -= satCharge * (*(shapedOvershootStart+index));
          *(shapedStart + sampleIdx) -= satAmp    * (*(shapedRingingStart+index));
        }
      }
    }
  } // End addition of pulses to raw signal

  // Convert raw signals to TDC, hit, ADC, Esum values
  // TDC values
  sample( m_rawSignal, m_tdcSignal,
           m_simFrequency, DayaBay::TdcFrequencyHz );
  std::vector<int> tdcValues;

  discriminate( m_tdcSignal, feeSimData->m_channelThreshold, tdcValues, 
        DayaBay::Threshold::kAbove);
  channel.setTdc(tdcValues);
  verbose() << "Channel " << channelId << ": " << tdcValues.size() 
        << " TDC values." << endreq;
  
  // Assemble Hit signal (for trigger) based on TDC clock values
  // The hit signal is a binary 0 or 1 vs. time at the Nhit frequency.
  unsigned int hitSamples = convertClock( m_tdcSignal.size(),
                      DayaBay::TdcFrequencyHz,
                      DayaBay::NhitFrequencyHz ) + 1;
  // Resize and clear the buffer
  if( m_hitSignal.size() != hitSamples) m_hitSignal.resize( hitSamples );
  int* hitStart = &m_hitSignal[0];
  for( unsigned int hitIdx = 0; hitIdx != hitSamples; hitIdx++ ){
    *(hitStart+hitIdx) = 0;
  }
  if( tdcValues.size() > 0 ){
    std::vector<int>::const_iterator tdcIter, tdcDone = tdcValues.end();
    for(tdcIter = tdcValues.begin(); tdcIter != tdcDone; tdcIter++){
      int tdcClock = *tdcIter;
      int hitClock = convertClock( tdcClock, DayaBay::TdcFrequencyHz, 
                   DayaBay::NhitFrequencyHz );
      m_hitSignal[hitClock] = 1;
    }
  }
  channel.setHit(m_hitSignal);
  verbose() << "Hit signal size: " << m_hitSignal.size() << endreq;
    
  // ADC Values: Digitized signals vs. time for both the high and low gains
  // Reduce shaped pulse to ADC frequency
  sample( m_shapedSignal, m_adcAnalogSignal,
      m_simFrequency, DayaBay::AdcFrequencyHz );
  // 12-bit ADC after shaped signals. One for adcHigh and one for adcLow.
  int adcBits = 12;
  digitize( m_adcAnalogSignal, m_adcHighGain,
        feeSimData->m_adcRangeHigh, adcBits, 
        feeSimData->m_adcBaselineHigh );
  digitize( m_adcAnalogSignal, m_adcLowGain,
        feeSimData->m_adcRangeLow, adcBits, 
        feeSimData->m_adcBaselineLow );
  channel.setAdcHigh( m_adcHighGain );
  channel.setAdcLow( m_adcLowGain );
  verbose() << "ADC signal size: " << m_adcHighGain.size() << endreq;
  int maxAdc = 0;
  for (unsigned int jj=0;jj<m_adcHighGain.size();jj++){
    if( m_adcHighGain[jj]>maxAdc) maxAdc = m_adcHighGain[jj];
  }
  // Channel energy: contribution to the Esum signal 
  channel.setEnergy(m_rawSignal);
  //channel.setEnergy(m_shapedSignal);
  return StatusCode::SUCCESS;
}  // End generate One channel
int EsIdealFeeTool::convertClock ( int  clock,
int  inputFrequency,
int  outputFrequency 
) [private]

Definition at line 1107 of file EsIdealFeeTool.cc.

                                                                                  {
  // Convert a clock cycle from one frequency to another.  Round down
  // for incoherent frequencies.
  return int( ( clock / double(inputFrequency) ) * outputFrequency );
}
DayaBay::ESumComp::ESumComp_t EsIdealFeeTool::boardEsumType ( int  boardId) [private]

Definition at line 882 of file EsIdealFeeTool.cc.

                                                                  {
  // FIXME: This should probably be moved into a service.
  //         Probably into the CableMappingSvc
  if (boardId >= 1 && boardId < 7)
    {
      return DayaBay::ESumComp::kESumLow;
    }

  if (boardId >= 7 && boardId < 13)
    {
      return DayaBay::ESumComp::kESumHigh;
    }

  if (boardId == 13)
    {
      return DayaBay::ESumComp::kESumNone;
    }
      
  return DayaBay::ESumComp::kUnknown;
}
void EsIdealFeeTool::convolve ( const DayaBay::AnalogSignal signal,
const DayaBay::AnalogSignal response,
DayaBay::AnalogSignal output 
) [private]

Definition at line 904 of file EsIdealFeeTool.cc.

{
  // Primitive convolution of signal and response, to produce
  // output.  Assumes values outside of range are zero.  A few simple
  // optimizations are use to speed up the process.
  int sN = signal.size();
  int rN = response.size();
  output.resize(sN+rN);
  const double* sigD = &(signal[0]);
  const double* resD = &(response[0]);
  double* outD = &(output[0]);
  for (int i=0; i<(sN+rN); i++) {
    double sum = 0;
    for (int j=0; j<rN; j++) {
      if (resD[j]==0)
        continue;
      if (i-j>=0 && i-j<sN)
        sum+=sigD[i-j]*resD[j];
    }
    outD[i]=sum;
  }
}
void EsIdealFeeTool::shape ( const DayaBay::AnalogSignal rawSignal,
DayaBay::AnalogSignal shapedSignal,
const DayaBay::AnalogSignal response,
std::map< int, int >  convolutionWindows 
) [private]

Definition at line 930 of file EsIdealFeeTool.cc.

                                                                                        {
  std::map<int,int>::iterator convIt = convolutionWindows.begin();
  int convStart, convEnd;
  shapedSignal.clear();
  int samples = rawSignal.size();
  shapedSignal.resize(samples, 0);
  DayaBay::AnalogSignal rawFrag;
  DayaBay::AnalogSignal shapedFrag;
  for(; convIt != convolutionWindows.end(); convIt++){
    convStart = convIt->first * 96;
    convEnd = convIt->second * 96;
    if (convEnd > samples) convEnd = samples;
    rawFrag.resize( (convEnd - convStart) );
    shapedFrag.resize( (convEnd - convStart) );
    std::copy( rawSignal.begin() + convStart, rawSignal.begin() + convEnd , rawFrag.begin() );
    convolve( rawFrag, response, shapedFrag );
    shapedFrag.resize( (convEnd - convStart) );
    std::copy( shapedFrag.begin(), shapedFrag.end(), shapedSignal.begin() + convStart );
    debug() << "Convolve from " << convStart *1.5625 << " ns to " << convEnd * 1.5625 << " ns" << endreq;
  }
}
void EsIdealFeeTool::digitize ( const DayaBay::AnalogSignal signal,
DayaBay::DigitalSignal output,
double  range,
int  bits,
double  offset 
) [private]

Definition at line 1012 of file EsIdealFeeTool.cc.

{
  // Simple digitization algorithm
  if(output.size() != signal.size()) output.resize(signal.size());
  int maxADC = int( pow(2, bits) );
  int signalN = signal.size();
  const double* inputStart = &signal[0];
  int* outputStart = &output[0];
  for (int i=0; i<signalN; i++) {
    int* curOutput = outputStart + i;
    *(curOutput) = int(((*(inputStart + i))/range)*maxADC + offset);
    if (*curOutput>maxADC)
      *curOutput=maxADC;
    if (*curOutput<0)
      *curOutput=0;
    /* Slower, but easier to read version of the same calculation
    outputStart[i] = int((signal[i]/range)*maxADC + offset);
    if (output[i]>maxADC)
      output[i]=maxADC;
    if (output[i]<0)
      output[i]=0;
    */
  }
}
void EsIdealFeeTool::discriminate ( const DayaBay::AnalogSignal signal,
double  threshold,
std::vector< int > &  output,
DayaBay::Threshold::Threshold_t  mode 
) [private]

Definition at line 1039 of file EsIdealFeeTool.cc.

{
  // Take input signal array and discriminate using threshold.  Return
  // vector of indices which satisfy the discriminator.
  //  Current Modes:
  //    - Threshold crossed
  //    - Above threshold
  threshold *= m_speAmp;
  threshold *= m_gainFactor;
  threshold *= m_discThreshScale;
  output.clear();
  const double* signalStart = &signal[0];
  unsigned int nSamples = signal.size();
  float maxS = 0;
  bool trigger=false;
  if (DayaBay::Threshold::kCrossing == mode) {
    bool aboveThreshold = false;
    for (unsigned int i=0;i < nSamples; ++i) {
      if (!aboveThreshold && (*(signalStart+i))>=threshold) {
        output.push_back(i);
        aboveThreshold = true;
      } else if (aboveThreshold && (*(signalStart+i))<threshold)
        aboveThreshold=false;
    }
  } else if (DayaBay::Threshold::kAbove == mode) {
    for (unsigned int i=0;i < nSamples; ++i){
      if((*(signalStart+i))>maxS) maxS= (*(signalStart+i));
      if ((*(signalStart+i))>=threshold){
        output.push_back(i);
        trigger = true;
      }
    }
  } else {
    error() << "discriminate(): Unknown Discrimination Mode" << endreq;
  }
}
double EsIdealFeeTool::pmtPulse ( double  deltaT,
int  nPulse 
) [private]

Definition at line 1132 of file EsIdealFeeTool.cc.

                                                         {
  // Return ideal single pe pulse with amplitude 1 V
  double width; // pulse width parameter
  double mu; // parameter determining the degree of asymmetry
  if ( nPulse > 1 ){
    width = getPulseWidth( nPulse); 
    mu    = getPulseForm( nPulse);
  }
  else {
    width = 7.5e-9;
    mu    = 0.45;
  }
  double shift = 6e-9 -width/1.5;
  if (deltaT-shift<0) return 0.;
  return - exp( -pow( log( (deltaT-shift)/width),2) 
               / ( 2*pow(mu,2) ) ) * Gaudi::Units::volt;
}
double EsIdealFeeTool::overshoot ( double  deltaT) [private]

Definition at line 1149 of file EsIdealFeeTool.cc.

                                              {
  if (deltaT < 0) return 0.;
  double amp = 0.045; // Relative overshoot amplitude for spe pulses
  // Fermi onset
  double t0   = 50e-9; 
  double t1   = 10e-9;
  double fermi = 1. / (exp( (t0 - deltaT) / t1) + 1.);
  // Exponential overshoot component
  double tau = 145.e-9; // Overshoot decay time in s
  double expoOS = exp(-(deltaT-87e-9)/tau);
  // Slower overshoot component
  double mean = 0.4e-6;
  double sigma = 0.08e-6;
  double t = deltaT -mean;
  double gausOS = 0.12 * exp(pow(t,2)/(-2*pow(sigma,2)));
  // Undershoot 
  mean = 0.65e-6;
  sigma = 0.12e-6;
  t = deltaT -mean;
  double undershoot = -0.03 * exp(pow(t,2)/(-2*pow(sigma,2)));
  
  return  amp * fermi * (expoOS + gausOS + undershoot)
          * Gaudi::Units::volt;
}
double EsIdealFeeTool::adcOvershoot ( double  deltaT) [private]

Definition at line 1174 of file EsIdealFeeTool.cc.

                                                 {
  if (deltaT < 0) return 0.;
  // Relative overshoot amplitude
  double amp = 0.015; 
  // Fermi onset
  double t0   = 130e-9; 
  double t1   =10e-9;
  double fermi = 0.7 / (exp( (t0 - deltaT) / t1) + 1.);
  // exponential
  double tau = 500.e-9; 
  double expo = 1.5*exp(-(deltaT)/tau);
  // add 2 Gaussian peaks
  double mean = 0.45e-6;
  double sigma = 0.12e-6;
  double gaus1 = 0.45 * exp(pow(deltaT-mean,2)/(-2*pow(sigma,2)));
  mean = 2e-6;
  sigma = 0.7e-6;
  double gaus2 = 0.12 * exp(pow(deltaT-mean,2)/(-2*pow(sigma,2)));

  return amp * (fermi * expo + gaus1 + gaus2)
         * Gaudi::Units::volt;
}
double EsIdealFeeTool::ringing ( double  deltaT) [private]

Definition at line 1197 of file EsIdealFeeTool.cc.

                                           {
  double phase = 1.1e-6; //phase shift in raw signal
  return ringingFast(deltaT) + ringingSlow(deltaT,phase);
}
double EsIdealFeeTool::adcRinging ( double  deltaT) [private]

Definition at line 1202 of file EsIdealFeeTool.cc.

                                              {
  double phase = 3.6e-6; //phase shift in shaped signal
  return ringingFast(deltaT) - ringingSlow(deltaT,phase);
}
double EsIdealFeeTool::ringingFast ( double  deltaT) [private]

Definition at line 1206 of file EsIdealFeeTool.cc.

                                               {
  double amp = 0.0016; // relative amplitude
  double tau = 0.8e-6; // decay time
  double t0 = 0.45e-6; // oscillation period
  double phase = 0.55e-6; 
  double t = deltaT - phase;
  if (t<0) return 0;
  return -amp * exp(-t/tau) * sin(2*3.14 *t/t0 )
          * Gaudi::Units::volt;
}
double EsIdealFeeTool::ringingSlow ( double  deltaT,
double  phase 
) [private]

Definition at line 1216 of file EsIdealFeeTool.cc.

                                                             {
  double amp = 0.0007; // relative amplitude
  double tau = 2.4e-6; // decay time
  double t0 = 3.6e-6;  // oscillation period
  double t = deltaT - phase;
  if (t<0) return 0;
  return amp * exp(-t/tau) * sin(2*3.14 *t/t0 )
         * Gaudi::Units::volt;
}
double EsIdealFeeTool::ampSatFactor ( int  nPulse) [private]

Definition at line 1113 of file EsIdealFeeTool.cc.

                                             {
  double q = double(nPulse);
  double qSat = 533.98;
  double a = 2.884;
  return saturationModel(q, qSat, a);
}
double EsIdealFeeTool::chargeSatFactor ( int  nPulse) [private]

Definition at line 1119 of file EsIdealFeeTool.cc.

                                                {
  double q = double(nPulse);
  double qSat = 939.77;
  double a = 2.113;
  return saturationModel(q, qSat, a);
}
double EsIdealFeeTool::saturationModel ( double  q,
double  qSat,
double  a 
) [private]

Definition at line 1125 of file EsIdealFeeTool.cc.

                                                                     {
  if(q<1) return 1;
  double num = pow(1 + 8*pow(q/qSat,a),0.5) - 1;
  double denom = 4* pow(q/qSat,a);
  return pow(num/denom,0.5);
}
double EsIdealFeeTool::pulseShaping ( double  deltaT) [private]

Definition at line 1079 of file EsIdealFeeTool.cc.

                                                 {
  // Gaussian CR-(RC)^4 pulse shaping
  // See Knoll, Radiation Detection and Measurement
  // Returns the pulse amplitude at time deltaT after a positive step
  // function of amplitude 1
  if (deltaT<0)
    return 0.;
  double t0 = 25.0e-9; // RC time constant in seconds
  double r = deltaT/t0;
  return pow( r, 4 ) * exp( -r );
}
double EsIdealFeeTool::esumShaping ( double  deltaT) [private]

Definition at line 1091 of file EsIdealFeeTool.cc.

                                                {
  //ESum shaping function
  if (deltaT<0)
    return 0.;
  double t0 = 56.0; // RC time constant in seconds
  double r = deltaT/t0;
  return exp( -r )*(1/t0);
}
double EsIdealFeeTool::getPulseWidth ( int  nPulse) [private]

Definition at line 1226 of file EsIdealFeeTool.cc.

                                              {
  // Parameter determining the width of the spe pulse
  double x  = double(nPulse);
  double p0 = 2.92682e+01;
  double p1 = 8.71668e+02;
  double p2 = 7.12771e+02;
  double p3 = 2.98393e-04;
  
  return (p0/(exp((p1-x)/p2)+1.)+p3*x)* 1e-9;
}
double EsIdealFeeTool::getPulseForm ( int  nPulse) [private]

Definition at line 1237 of file EsIdealFeeTool.cc.

                                             {
  // Parameter determining the asymmetry of the spe pulse
  double x = double(nPulse);
  double p0 =  2.03759e-01;
  double p1 = -8.07750e+02;
  double p2 =  5.18878e-02;
  double p3 = -7.48266e+02;
  double p4 =  2.19051e-01;
  
  return p0*exp(x/p1) + p2*exp(x/p3) + p4;
}
double EsIdealFeeTool::getPulseShift ( int  nPulse) [private]
double EsIdealFeeTool::getOvershootAmp ( int  nPulse) [private]
double EsIdealFeeTool::noise ( ) [private]

Definition at line 1100 of file EsIdealFeeTool.cc.

                             {
  // FIXME: do discrete FT for range of frequencies
  // Simple noise model
  // Approximate Gaussian-distributed noise
  return m_gauss() * m_noiseAmp; // no offset
}
void EsIdealFeeTool::sample ( const DayaBay::AnalogSignal signal,
DayaBay::AnalogSignal output,
int  inputFrequency,
int  outputFrequency 
) [private]

Definition at line 953 of file EsIdealFeeTool.cc.

{
  // Fill the output array using the signal values sampled at a
  // regular frequency
  if(inputFrequency == outputFrequency){
    if( output.size() != signal.size() ) output.resize( signal.size() );
    memcpy(&output[0],&signal[0],output.size()*sizeof(double));
    return;
  }
  unsigned int nSamples = convertClock(signal.size(), 
                       inputFrequency, 
                       outputFrequency);
  if(output.size() != nSamples){
    if(output.capacity() > nSamples){
      output.resize(nSamples);
    }else{
      // Faster memory allocation
      DayaBay::AnalogSignal* newSignal = new DayaBay::AnalogSignal(nSamples);
      output.swap(*newSignal);
      delete newSignal;
    }
  }

  if(inputFrequency < outputFrequency){
    warning() << "sample(): Input frequency " << inputFrequency 
          << " is less than the sampling frequency " << outputFrequency
          << endreq;
  }
  if( (inputFrequency % outputFrequency) == 0 ){
    // Check for integer relationship between frequencies
    int relativeFrequency = inputFrequency / outputFrequency;
    double* outputStart = &output[0]; 
    const double* inputStart = &signal[0]; 
    int inputIdx = 0;
    unsigned int outputIdx = 0;
    while(outputIdx != nSamples){
      *(outputStart + outputIdx) = *(inputStart + inputIdx);
      outputIdx++;
      inputIdx+=relativeFrequency;
    }
    return;
  }else{
    // Frequencies not coherent.  Use signal sample with nearest lower index.
    double* outputStart = &output[0];
    const double* inputStart = &signal[0];
    int inputIdx = 0;
    unsigned int outputIdx = 0;
    double freqScale = inputFrequency / outputFrequency;
    while( outputIdx != nSamples ){
      inputIdx = int(outputIdx * freqScale);
      *(outputStart + outputIdx) = *(inputStart + inputIdx);
      outputIdx++;
    }
  }
}
const InterfaceID & IEsFrontEndTool::interfaceID ( ) [static, inherited]

Retrieve interface ID.

Definition at line 8 of file IEsFrontEndTool.cc.

{ 
    return IID_IEsFrontEndTool; 
}

Member Data Documentation

std::string EsIdealFeeTool::m_cableSvcName [private]

Definition at line 108 of file EsIdealFeeTool.h.

std::string EsIdealFeeTool::m_simDataSvcName [private]

Definition at line 110 of file EsIdealFeeTool.h.

Definition at line 112 of file EsIdealFeeTool.h.

Definition at line 114 of file EsIdealFeeTool.h.

Definition at line 116 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_gainFactor [private]

Definition at line 118 of file EsIdealFeeTool.h.

Definition at line 120 of file EsIdealFeeTool.h.

Definition at line 122 of file EsIdealFeeTool.h.

Definition at line 124 of file EsIdealFeeTool.h.

unsigned int EsIdealFeeTool::m_linearityThreshold [private]

Definition at line 126 of file EsIdealFeeTool.h.

Definition at line 129 of file EsIdealFeeTool.h.

Definition at line 131 of file EsIdealFeeTool.h.

Definition at line 132 of file EsIdealFeeTool.h.

Definition at line 134 of file EsIdealFeeTool.h.

Definition at line 136 of file EsIdealFeeTool.h.

Definition at line 138 of file EsIdealFeeTool.h.

Definition at line 140 of file EsIdealFeeTool.h.

Definition at line 142 of file EsIdealFeeTool.h.

Definition at line 144 of file EsIdealFeeTool.h.

Definition at line 146 of file EsIdealFeeTool.h.

Definition at line 147 of file EsIdealFeeTool.h.

Definition at line 148 of file EsIdealFeeTool.h.

Definition at line 150 of file EsIdealFeeTool.h.

Definition at line 152 of file EsIdealFeeTool.h.

Definition at line 154 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_noiseAmp [private]

Definition at line 156 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_speAmp [private]

Definition at line 158 of file EsIdealFeeTool.h.

Definition at line 160 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_adcRange [private]

Definition at line 162 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_adcBits [private]

Definition at line 163 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_adcOffset [private]

Definition at line 164 of file EsIdealFeeTool.h.

Definition at line 166 of file EsIdealFeeTool.h.

Definition at line 168 of file EsIdealFeeTool.h.

Definition at line 170 of file EsIdealFeeTool.h.

Definition at line 172 of file EsIdealFeeTool.h.

Definition at line 174 of file EsIdealFeeTool.h.

Definition at line 176 of file EsIdealFeeTool.h.

Definition at line 178 of file EsIdealFeeTool.h.

Definition at line 180 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_shapingSum [private]

Definition at line 182 of file EsIdealFeeTool.h.

Definition at line 184 of file EsIdealFeeTool.h.

double EsIdealFeeTool::m_pulseMax [private]

Definition at line 186 of file EsIdealFeeTool.h.

Rndm::Numbers EsIdealFeeTool::m_gauss [private]

Definition at line 189 of file EsIdealFeeTool.h.

Definition at line 191 of file EsIdealFeeTool.h.

Definition at line 192 of file EsIdealFeeTool.h.

Definition at line 193 of file EsIdealFeeTool.h.

Definition at line 194 of file EsIdealFeeTool.h.

Definition at line 195 of file EsIdealFeeTool.h.

Definition at line 196 of file EsIdealFeeTool.h.

Definition at line 197 of file EsIdealFeeTool.h.

Definition at line 198 of file EsIdealFeeTool.h.

Definition at line 199 of file EsIdealFeeTool.h.

Definition at line 200 of file EsIdealFeeTool.h.

Definition at line 201 of file EsIdealFeeTool.h.

Definition at line 202 of file EsIdealFeeTool.h.

Definition at line 203 of file EsIdealFeeTool.h.

Definition at line 204 of file EsIdealFeeTool.h.

Definition at line 205 of file EsIdealFeeTool.h.

Definition at line 206 of file EsIdealFeeTool.h.

Definition at line 207 of file EsIdealFeeTool.h.

Definition at line 208 of file EsIdealFeeTool.h.

std::map<int,int> EsIdealFeeTool::m_adcConvolutionWindows [private]

Definition at line 209 of file EsIdealFeeTool.h.

std::map<int,int> EsIdealFeeTool::m_esumConvolutionWindows [private]

Definition at line 210 of file EsIdealFeeTool.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:16:37 for ElecSim by doxygen 1.7.4