00001 #include "EsIdealFeeTool.h"
00002 #include "Event/ElecFeeCrate.h"
00003 #include "Event/ElecPulseCollection.h"
00004 #include "Event/ElecPmtPulse.h"
00005 #include "GaudiKernel/SystemOfUnits.h"
00006 #include "Conventions/Trigger.h"
00007
00008 #include <sstream>
00009 #include <string.h>
00010
00011 EsIdealFeeTool::EsIdealFeeTool(const std::string& type,
00012 const std::string& name,
00013 const IInterface* parent)
00014 : GaudiTool(type,name,parent),
00015 m_cableSvc(0)
00016 {
00017 declareInterface< IEsFrontEndTool >(this);
00018
00019 declareProperty("CableSvcName",m_cableSvcName="StaticCableSvc",
00020 "Name of service to map between detector, hardware, and electronic IDs");
00021 declareProperty("SimDataSvcName",m_simDataSvcName="StaticSimDataSvc",
00022 "Name of service to provide FEE channel properties for simulation");
00023 declareProperty("SimFrequency", m_simFrequency = DayaBay::TdcFrequencyHz,
00024 "Simulation Frequency");
00025 declareProperty("ESumFrequency", m_eSumFrequency = DayaBay::EsumFrequencyHz,
00026 "Esum Frequency");
00027 declareProperty("TriggerWindowCycles",
00028 m_triggerWindowCycles=DayaBay::TriggerWindowCycles,
00029 "Number of trigger window ccyles");
00030 declareProperty("GainFactor",
00031 m_gainFactor=2.,
00032 "Mean PMT gain (in units of 10^7)");
00033 declareProperty("NoiseAmp",m_noiseAmp=0.5e-3*Gaudi::Units::volt,
00034 "Voltage amplitude of noise");
00035 declareProperty("SpeAmp",m_speAmp=0.0035,
00036 "Spe pulse height in V at 1e7 gain");
00037 declareProperty("DiscThreshScale", m_discThreshScale=1,
00038 "Factor to scale all channels' discriminator threshold");
00039 declareProperty("EnableNonlinearity", m_enableNonlinearity=true,
00040 "Turns on/off all nonlinear effects");
00041 declareProperty("EnableNoise",m_enableNoise=false,"Turn on/off noise individually");
00042 declareProperty("EnableDynamicWaveform", m_enableDynamicWaveform=true,
00043 "Turns on/off dynamic waveform individually");
00044 declareProperty("EnableRinging", m_enableRinging=true,
00045 "Turns on/off ringing individually");
00046 declareProperty("EnableOvershoot", m_enableOvershoot=true,
00047 "Turns on/off overshoot individually");
00048 declareProperty("EnableSaturation", m_enableSaturation=true,
00049 "Turns on/off saturation individually");
00050 declareProperty("EnableAcCoupling", m_enableAcCoupling=true,
00051 "Turns on/off high pass filter before discriminator");
00052 declareProperty("EnableESumTotal", m_enableESumTotal=true,
00053 "Turns on/off Total ESum simulation");
00054 declareProperty("EnableESumH", m_enableESumH=true, "Turns on/off upper ESum simulation");
00055 declareProperty("EnableESumL", m_enableESumL=true, "Turns on/off lower ESum simulation");
00056 declareProperty("EnableFastSimMode", m_enableFastSimMode=false,
00057 "Turns on/off fast simulation mode to save CPU time");
00058 declareProperty("PulseCountSlotWidth", m_pulseCountSlotWidth=5e-9*Gaudi::Units::second,
00059 "pulse count time slot size");
00060 declareProperty("PulseCountWindow", m_pulseCountWindow=2,
00061 "Number of counted slots before and after a pulse");
00062 declareProperty("LinearityThreshold", m_linearityThreshold=20,
00063 "Threshold in P.E. for activating nonlinear pulse shape");
00064 declareProperty("EnablePretrigger", m_enablePretrigger=true,
00065 "Turns on/off pretrigger");
00066 declareProperty("TotalESumTriggerThreshold",
00067 m_eSumTriggerThreshold=DayaBay::Trigger::kADESumThreshold,
00068 "Total ESum trigger threshold as set in TrigSim");
00069 declareProperty("NHitTriggerThreshold",
00070 m_nHitTriggerThreshold=DayaBay::Trigger::kADthreshold,
00071 "nHit trigger threshold as set in TrigSim");
00072 declareProperty("ADCOffset", m_adcOffset = 500, "Esum Frequency");
00073 declareProperty("ADCRange", m_adcRange = 1, "Esum Range");
00074 declareProperty("ADCBits", m_adcBits = 10, "Esum Bits");
00075 }
00076
00077 EsIdealFeeTool::~EsIdealFeeTool(){}
00078
00079 StatusCode EsIdealFeeTool::initialize()
00080 {
00081
00082 m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
00083
00084
00085 m_simDataSvc = svc<ISimDataSvc>(m_simDataSvcName,true);
00086
00087 StatusCode sc = loadResponse();
00088 if(sc.isFailure()) return sc;
00089
00090
00091 IRndmGenSvc *rgs = 0;
00092 if (service("RndmGenSvc",rgs,true).isFailure()) {
00093 fatal() << "Failed to get random service" << endreq;
00094 return StatusCode::FAILURE;
00095 }
00096
00097 sc = m_gauss.initialize(rgs, Rndm::Gauss(0,1));
00098 if (sc.isFailure()) {
00099 fatal() << "Failed to initialize gaussian random numbers" << endreq;
00100 return StatusCode::FAILURE;
00101 }
00102
00103 int preThr = int(0.8 * m_nHitTriggerThreshold + 0.5);
00104 int esumPreThr = int (0.8 * 1500 * m_eSumTriggerThreshold/Gaudi::Units::volt + 0.5);
00105 if(preThr > esumPreThr) preThr = esumPreThr;
00106 m_pretriggerThreshold = preThr;
00107 if( m_enablePretrigger ){
00108 info() << "Pretrigger threshold is " << m_pretriggerThreshold << " pe " << endreq;
00109 }
00110 else {
00111 info() << "Pretrigger disabled" << endreq;
00112 info() << "Convolution threshold is " << m_pretriggerThreshold << " pe " << endreq;
00113 }
00114
00115 if( m_enableFastSimMode){
00116 info() << "Fast electronics simulation mode enabled: Note that Upper and Lower ESum trigger do not yield meaningful results" << endreq;
00117 m_linearityThreshold = 100000;
00118 m_enableESumH = false;
00119 m_enableESumL = false;
00120 }
00121 if( !m_enableNonlinearity ) {
00122 info() << "All nonlinear effects disabled." << endreq;
00123 m_enableOvershoot = false;
00124 m_enableSaturation = false;
00125 m_enableRinging = false;
00126 m_enableDynamicWaveform = false;
00127 m_enableAcCoupling = false;
00128 }
00129 else {
00130 info() << "Nonlinear model enabled for channels with more than "
00131 << m_linearityThreshold << " hits " << endreq;
00132
00133 if( !m_enableAcCoupling )
00134 info() << "AC coupling disabled" << endreq;
00135 else
00136 debug() << "AC coupling enabled" << endreq;
00137
00138 if( !m_enableOvershoot )
00139 info() << "Overshoot disabled" << endreq;
00140 else
00141 debug() << "Overshoot enabled" << endreq;
00142
00143 if( !m_enableRinging )
00144 info() << "Ringing disabled" << endreq;
00145 else
00146 debug() << "Ringing enabled" << endreq;
00147
00148 if( !m_enableSaturation )
00149 info() << "Saturation disabled" << endreq;
00150 else
00151 debug() << "Saturation enabled" << endreq;
00152
00153 if( !m_enableDynamicWaveform )
00154 info() << "Dynamic waveform disabled" << endreq;
00155 else
00156 debug() << "Dynamic waveform enabled" << endreq;
00157 }
00158 return StatusCode::SUCCESS;
00159 }
00160
00161 StatusCode EsIdealFeeTool::finalize()
00162 {
00163 return StatusCode::SUCCESS;
00164 }
00165
00166 StatusCode EsIdealFeeTool::generateSignals(DayaBay::ElecPulseCollection* pulses,
00167 DayaBay::ElecCrate* elecCrate)
00168 {
00169
00170 DayaBay::ElecFeeCrate* crate = 0;
00171 crate = dynamic_cast<DayaBay::ElecFeeCrate*>(elecCrate);
00172 if(crate == NULL){
00173 error() << "generateSignals(): Crate is not a FEE Crate." << endreq;
00174 return StatusCode::FAILURE;
00175 }
00176
00177
00178 Context context(pulses->detector().site(), SimFlag::kMC,
00179 pulses->header()->header()->timeStamp(),
00180 pulses->detector().detectorId());
00181 int task = 0;
00182 ServiceMode svcMode(context, task);
00183
00184
00185 if(crate->channelData().size() == 0){
00186
00187 const std::vector<DayaBay::DetectorSensor>& pmtList
00188 = m_cableSvc->sensors( svcMode );
00189 int nPmts = pmtList.size();
00190 for(int pmtIdx = 0; pmtIdx < nPmts; pmtIdx++){
00191
00192 DayaBay::FeeChannelId channelId(m_cableSvc->elecChannelId(pmtList[pmtIdx], svcMode).fullPackedData());
00193 crate->addChannel(channelId);
00194 verbose() << "Adding channel to crate: " << channelId << endreq;
00195 }
00196 }
00197
00198
00199 TimeStamp earliest = pulses->header()->header()->earliest();
00200 TimeStamp latest = pulses->header()->header()->latest();
00201
00202 double simTime = latest - earliest;
00203 verbose() << "Simulation time is " << simTime << endreq;
00204
00205 int simSamples = int(simTime * m_simFrequency);
00206 verbose() << "Simulation has " << simSamples << " samples at frequency "
00207 << m_simFrequency << " Hz" << endreq;
00208
00209
00210 int nTimeSlices = int(simTime / 150e-9 + 0.5);
00211
00212
00213 int hitCountSlices[nTimeSlices];
00214 for ( int i = 0; i < nTimeSlices; ++i ) hitCountSlices[i] = 0;
00215
00216
00217 std::map<DayaBay::ElecChannelId, std::vector<DayaBay::ElecPmtPulse*> > pulseMap;
00218 StatusCode sc;
00219 sc = mapPulsesByChannel( pulses->pulses(), pulseMap, hitCountSlices );
00220 if(sc.isFailure()) return sc;
00221 debug() << "Mapped Pulses" << endreq;
00222
00223
00224
00225 m_adcConvolutionWindows.clear();
00226 m_esumConvolutionWindows.clear();
00227 int adcConvStart, adcConvEnd;
00228 int esumConvStart, esumConvEnd;
00229 std::map<int,int>::reverse_iterator sliceIt;
00230 for(int i = 1; i < nTimeSlices; i++){
00231 if( hitCountSlices[i] + hitCountSlices[i-1] < m_pretriggerThreshold ) continue;
00232 if (i - 5 < 0) adcConvStart = 0;
00233 else adcConvStart = (i - 5);
00234 if (i - 1 < 0) esumConvStart = 0;
00235 else esumConvStart = (i - 1);
00236
00237 if (i + 5 > nTimeSlices) adcConvEnd = nTimeSlices;
00238 else adcConvEnd = (i + 5);
00239 if (i + 1 > nTimeSlices) esumConvEnd = nTimeSlices;
00240 else esumConvEnd = (i + 1);
00241
00242 if( m_adcConvolutionWindows.empty() ) {
00243 m_adcConvolutionWindows[adcConvStart] = adcConvEnd;
00244 m_esumConvolutionWindows[esumConvStart] = esumConvEnd;
00245 continue;
00246 }
00247 sliceIt = m_adcConvolutionWindows.rbegin();
00248 if( (adcConvStart - sliceIt->second > 2) )
00249 m_adcConvolutionWindows[adcConvStart ] = adcConvEnd;
00250 else
00251 m_adcConvolutionWindows[sliceIt->first] = adcConvEnd;
00252
00253 sliceIt = m_esumConvolutionWindows.rbegin();
00254 if( (esumConvStart - sliceIt->second > 2) )
00255 m_esumConvolutionWindows[esumConvStart ] = esumConvEnd;
00256 else
00257 m_esumConvolutionWindows[sliceIt->first] = esumConvEnd;
00258 }
00259 debug() << "Number of time windows to be possibly read out: " << m_adcConvolutionWindows.size() << endreq;
00260
00261
00262 if(m_adcConvolutionWindows.size()==0 && m_enablePretrigger) {
00263 debug() << "Event below pretrigger threshold, no waveform simulation" << endreq;
00264 return StatusCode::SUCCESS;
00265 }
00266
00267
00268 const DayaBay::ElecFeeCrate::ChannelData& channelData = crate->channelData();
00269 DayaBay::ElecFeeCrate::ChannelData::const_iterator channelIter,
00270 done = channelData.end();
00271 for(channelIter = channelData.begin(); channelIter != done; ++channelIter){
00272 DayaBay::FeeChannelId channelId = channelIter->first;
00273 DayaBay::ElecFeeChannel& channel = crate->channel(channelId);
00274
00275
00276 sc = generateOneChannel(channelId, pulseMap[channelId], channel,
00277 simTime, svcMode);
00278 if( !sc.isSuccess() ) return sc;
00279
00280 }
00281
00282 StatusCode status = updateBoardSignals(crate);
00283 if (status.isFailure()) return status;
00284 status = updateEsumSignals(crate);
00285 if (status.isFailure()) return status;
00286
00287 verbose() << "Crate\n" << *crate << endreq;
00288
00289 return StatusCode::SUCCESS;
00290 }
00291
00292 StatusCode EsIdealFeeTool::generateOneChannel(
00293 const DayaBay::FeeChannelId& channelId,
00294 std::vector<DayaBay::ElecPmtPulse*>& channelPulses,
00295 DayaBay::ElecFeeChannel& channel,
00296 double simTime,
00297 const ServiceMode& svcMode)
00298 {
00299
00300
00301
00302 const DayaBay::FeeSimData* feeSimData =
00303 m_simDataSvc->feeSimData(channelId, svcMode);
00304 if(!feeSimData){
00305 error() << "No Simulation input properties for FEE channel: "
00306 << channelId << endreq;
00307 return StatusCode::FAILURE;
00308 }
00309
00310
00311 unsigned int simSamples = int(simTime * m_simFrequency);
00312
00313
00314 if(m_rawSignal.size() != simSamples) m_rawSignal.resize( simSamples );
00315 if(m_discriminatorSignal.size() != simSamples) m_discriminatorSignal.resize( simSamples );
00316 if(m_shapedSignal.size() != simSamples) m_shapedSignal.resize( simSamples );
00317 double* rawStart = &m_rawSignal[0];
00318 double* discriminatorStart = &m_discriminatorSignal[0];
00319 double* shapedStart = &m_shapedSignal[0];
00320 for( unsigned int sigIdx = 0; sigIdx!=simSamples; sigIdx++){
00321 *(rawStart + sigIdx) = 0;
00322 *(discriminatorStart + sigIdx) = 0;
00323 *(shapedStart + sigIdx) = 0;
00324 }
00325
00326
00327 if( m_enableNoise ){
00328 for(unsigned int index=0; index < simSamples; index++){
00329 *(rawStart + index) += noise();
00330 *(discriminatorStart + index) += noise();
00331 }
00332 }
00333
00334
00335 if( channelPulses.size() > 0 ){
00336 verbose() << "Channel " << channelId << " has " << channelPulses.size()
00337 << " pulses." << endreq;
00338
00339 if( m_enableNonlinearity && (channelPulses.size() > m_linearityThreshold) ){
00340
00341
00342 verbose() << "Applying nonlinearity:" << channelPulses.size()
00343 << " pulses (threshold=" << m_linearityThreshold << ")"
00344 << endreq;
00345
00346
00347 int numPulseTimeSlots = int(simTime*Gaudi::Units::second
00348 /m_pulseCountSlotWidth) + 1;
00349 verbose() << "Number of time slots for afterpulse counting = "
00350 << numPulseTimeSlots << std::endl;
00351 std::vector<int> pulseTimeSlots(numPulseTimeSlots);
00352
00353
00354
00355 for (int i = 0; i < numPulseTimeSlots; i++ ) pulseTimeSlots[i] = 0;
00356 std::vector<DayaBay::ElecPmtPulse*>::const_iterator pulseIter, pulseDone
00357 = channelPulses.end();
00358 for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
00359 DayaBay::ElecPmtPulse* pulse = *pulseIter;
00360 int timeSlot = int(pulse->time() * Gaudi::Units::nanosecond
00361 / m_pulseCountSlotWidth);
00362 int pulseNo = int(pulse->amplitude()+0.5);
00363 pulseTimeSlots[timeSlot]+=pulseNo;
00364 verbose() << "Added " << pulseNo << " pulses in time slot " << timeSlot << endreq;
00365 }
00366
00367
00368 for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
00369 DayaBay::ElecPmtPulse* pulse = *pulseIter;
00370 int tOffset = int(pulse->time()*1e-9 * m_simFrequency);
00371 float amplitude = pulse->amplitude();
00372
00373
00374
00375
00376 int nPulse=0;
00377 int timeSlot = int(pulse->time() * Gaudi::Units::nanosecond
00378 / m_pulseCountSlotWidth);
00379 for (int i = timeSlot - m_pulseCountWindow;
00380 i <= timeSlot + m_pulseCountWindow; i++){
00381 if(i>=0 && i<numPulseTimeSlots) nPulse += pulseTimeSlots[i];
00382 }
00383 verbose() << "Number of pulses in time window: " << nPulse << endreq;
00384
00385
00386 int nPulseSamples = m_pmtPulse.size();
00387 for(int index = 0; index < nPulseSamples; index++){
00388 unsigned int sampleIdx = tOffset + index;
00389 double idxTime = index * (1. / m_simFrequency);
00390 if(sampleIdx>0 && sampleIdx<simSamples){
00391
00392 if( m_enableDynamicWaveform ){
00393 m_rawSignal[sampleIdx] += - amplitude * pmtPulse( idxTime, nPulse, m_enableOvershoot );
00394 m_discriminatorSignal[sampleIdx] += - amplitude * pmtPulse( idxTime, nPulse, false );
00395 }else{
00396 if(m_enableOvershoot) m_rawSignal[sampleIdx] += - amplitude * m_pmtPulse[index];
00397 else m_rawSignal[sampleIdx] += - amplitude * m_pmtPulse_noOvershoot[index];
00398 m_discriminatorSignal[sampleIdx] += - amplitude * m_pmtPulse_noOvershoot[index];
00399 }
00400 }
00401 }
00402 }
00403
00404
00405 if( (channelPulses.size() > 0)
00406 && (m_enableSaturation || m_enableRinging) ){
00407
00408
00409
00410 double rawMax;
00411 int maxSample;
00412 rawMax = 0;
00413 for (unsigned int index=0; index<simSamples; index++) {
00414
00415 if( m_rawSignal[index] > rawMax ) {
00416 rawMax = m_rawSignal[index];
00417 maxSample = index;
00418 }
00419 if(m_enableSaturation){
00420 m_rawSignal[index] *= satFactor(m_rawSignal[index]);
00421 m_discriminatorSignal[index] *= satFactor(m_discriminatorSignal[index]);
00422 }
00423 }
00424
00425
00426 if( m_enableRinging && rawMax != 0 ){
00427 for (unsigned int index=maxSample; index<simSamples; index++) {
00428 double idxTime = index * (1. / m_simFrequency);
00429 m_rawSignal[index] += - ringing ( idxTime, rawMax );
00430 }
00431 }
00432 }
00433
00434
00435 shape(m_rawSignal,m_shapedSignal,m_shapingResponse,m_adcConvolutionWindows);
00436 m_shapedSignal.resize( simSamples );
00437
00438 for(unsigned int index=0; index < simSamples; index++) {
00439 m_shapedSignal[index] /= m_shapingSum;
00440 }
00441
00442
00443 }else{
00444
00445
00446 double* pmtPulseStart;
00447 if(m_enableOvershoot) pmtPulseStart = &m_pmtPulse[0];
00448 else pmtPulseStart = &m_pmtPulse_noOvershoot[0];
00449 double* pmtDiscriminatorPulseStart = &m_pmtPulse_noOvershoot[0];
00450 double* shapedPmtPulseStart = &m_shapedPmtPulse[0];
00451 double* ringingStart = &m_ringing[0];
00452 double* shapedRingingStart = &m_shapedRinging[0];
00453
00454
00455 std::vector<DayaBay::ElecPmtPulse*>::const_iterator pulseIter, pulseDone
00456 = channelPulses.end();
00457 for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
00458 DayaBay::ElecPmtPulse* pulse = *pulseIter;
00459 int tOffset = int(pulse->time()*1e-9 * m_simFrequency);
00460 float amplitude = pulse->amplitude();
00461
00462
00463 int nPulseSamples = m_pmtPulse.size();
00464 for(int index = 0; index < nPulseSamples; index++){
00465 unsigned int sampleIdx = tOffset + index;
00466
00467 if(sampleIdx>0 && sampleIdx<simSamples){
00468 *(rawStart + sampleIdx) -= amplitude * (*(pmtPulseStart + index));
00469 *(discriminatorStart + sampleIdx) -= amplitude * (*(pmtDiscriminatorPulseStart + index));
00470 *(shapedStart + sampleIdx) -= amplitude * (*(shapedPmtPulseStart+index));
00471 }
00472 }
00473 }
00474
00475
00476 if( (channelPulses.size() > 0) && (m_enableSaturation || m_enableRinging) ){
00477
00478
00479 double rawMax;
00480 int maxSample;
00481 rawMax = 0;
00482 for (unsigned int index=0; index<simSamples; index++) {
00483 double* curValue = rawStart + index;
00484 if( *curValue > rawMax ) {
00485 rawMax = *curValue;
00486 maxSample = index;
00487 }
00488 if(m_enableFastSimMode && m_enableSaturation){
00489 *(rawStart + index) *= satFactor(*(rawStart + index));
00490 *(rawStart + index) *= satFactor(-10 * (*(rawStart + index)));
00491 }
00492 }
00493
00494 if( m_enableRinging && rawMax != 0 ){
00495 unsigned int ringEnd = m_ringing.size();
00496 double ringScaling = -rawMax / m_pulseMax;
00497 if( ringEnd + maxSample > simSamples )
00498 ringEnd = simSamples - maxSample;
00499 for (unsigned int ringIdx=0; ringIdx<ringEnd; ringIdx++) {
00500 int index = maxSample + ringIdx;
00501 if(m_enableFastSimMode && m_enableSaturation){
00502 *(discriminatorStart + index) *= satFactor(*(discriminatorStart + index));
00503 *(shapedStart + index) *= pow(satFactor(15.* (fabs(*(shapedStart + index))) ),0.75);
00504 }
00505 *(rawStart + index) += ringScaling * (*(ringingStart+ringIdx));
00506 *(shapedStart + index) += ringScaling * (*(shapedRingingStart+ringIdx));
00507 }
00508 }
00509 }
00510
00511 }
00512
00513 }
00514
00515
00516
00517 if (m_enableAcCoupling) sample( m_discriminatorSignal, m_tdcSignal,
00518 m_simFrequency, DayaBay::TdcFrequencyHz );
00519 else sample( m_rawSignal, m_tdcSignal,
00520 m_simFrequency, DayaBay::TdcFrequencyHz );
00521 std::vector<int> tdcValues;
00522 discriminate( m_tdcSignal, feeSimData->m_channelThreshold, tdcValues,
00523 DayaBay::Threshold::kAbove);
00524 channel.setTdc(tdcValues);
00525 verbose() << "Channel " << channelId << ": " << tdcValues.size()
00526 << " TDC values." << endreq;
00527
00528
00529
00530 unsigned int hitSamples = convertClock( m_tdcSignal.size(),
00531 DayaBay::TdcFrequencyHz,
00532 DayaBay::NhitFrequencyHz ) + 1;
00533
00534 if( m_hitSignal.size() != hitSamples) m_hitSignal.resize( hitSamples );
00535 int* hitStart = &m_hitSignal[0];
00536 for( unsigned int hitIdx = 0; hitIdx != hitSamples; hitIdx++ ){
00537 *(hitStart+hitIdx) = 0;
00538 }
00539 if( tdcValues.size() > 0 ){
00540 std::vector<int>::const_iterator tdcIter, tdcDone = tdcValues.end();
00541 for(tdcIter = tdcValues.begin(); tdcIter != tdcDone; tdcIter++){
00542 int tdcClock = *tdcIter;
00543 int hitClock = convertClock( tdcClock, DayaBay::TdcFrequencyHz,
00544 DayaBay::NhitFrequencyHz );
00545 m_hitSignal[hitClock] = 1;
00546 }
00547 }
00548 channel.setHit(m_hitSignal);
00549 verbose() << "Hit signal size: " << m_hitSignal.size() << endreq;
00550
00551
00552
00553 sample( m_shapedSignal, m_adcAnalogSignal,
00554 m_simFrequency, DayaBay::AdcFrequencyHz );
00555
00556 int adcBits = 12;
00557 digitize( m_adcAnalogSignal, m_adcHighGain,
00558 feeSimData->m_adcRangeHigh, adcBits,
00559 feeSimData->m_adcBaselineHigh );
00560 digitize( m_adcAnalogSignal, m_adcLowGain,
00561 feeSimData->m_adcRangeLow, adcBits,
00562 feeSimData->m_adcBaselineLow );
00563 channel.setAdcHigh( m_adcHighGain );
00564 channel.setAdcLow( m_adcLowGain );
00565 verbose() << "ADC signal size: " << m_adcHighGain.size() << endreq;
00566
00567
00568
00569
00570
00571 channel.setEnergy(m_rawSignal);
00572 return StatusCode::SUCCESS;
00573 }
00574
00576
00577
00578 StatusCode EsIdealFeeTool::loadResponse(){
00579
00580 double dT_seconds = (1. / m_simFrequency);
00581 double dT_units = dT_seconds*CLHEP::second;
00582 int nShapingSamples = int(DayaBay::preTimeTolerance/dT_units);
00583 int nPulseSamples = 3*nShapingSamples;
00584 int nRingingSamples = int(DayaBay::postTimeTolerance/dT_units);
00585 m_shapingResponse.resize(nShapingSamples);
00586 m_esumResponse.resize(nPulseSamples);
00587 m_pmtPulse.resize(nPulseSamples);
00588 m_pmtPulse_noOvershoot.resize(nPulseSamples);
00589 m_ringing.resize(nRingingSamples);
00590 m_shapingSum = 0;
00591 m_esumShapingSum = 0;
00592 m_pulseMax = 0;
00593
00594 for (int i=0; i<nShapingSamples; i++) {
00595 m_shapingResponse[i] = pulseShaping(i*dT_seconds);
00596 m_shapingSum += m_shapingResponse[i];
00597 }
00598
00599 for (int i=0; i<nPulseSamples; i++){
00600 m_esumResponse[i] = esumShaping(i*dT_units/CLHEP::nanosecond);
00601 m_esumShapingSum += m_esumResponse[i];
00602 }
00603
00604 for (int i=0; i<nPulseSamples; i++) {
00605 m_pmtPulse[i] = pmtPulse(i*dT_seconds,1,true);
00606 m_pmtPulse_noOvershoot[i] = pmtPulse(i*dT_seconds,1,false);
00607
00608 if( m_pmtPulse[i] < m_pulseMax ) m_pulseMax = m_pmtPulse[i];
00609
00610 }
00611
00612 m_pulseMax *= -1;
00613
00614 for (int i=0; i<nRingingSamples; i++) {
00615 m_ringing[i] = ringing(i*dT_seconds, m_pulseMax);
00616 }
00617
00618
00619 convolve( m_pmtPulse, m_shapingResponse, m_shapedPmtPulse );
00620 m_shapedPmtPulse.resize( m_pmtPulse.size() );
00621 convolve( m_pmtPulse_noOvershoot, m_shapingResponse, m_shapedPmtPulse_noOvershoot );
00622 m_shapedPmtPulse_noOvershoot.resize( m_pmtPulse.size() );
00623
00624 for(unsigned int i=0; i < m_pmtPulse.size(); i++) {
00625 m_shapedPmtPulse[i] /= m_shapingSum;
00626 m_shapedPmtPulse_noOvershoot[i] /= m_shapingSum;
00627 }
00628
00629 convolve( m_ringing, m_shapingResponse, m_shapedRinging );
00630 m_shapedRinging.resize( m_ringing.size() );
00631
00632 for(unsigned int i=0; i < m_shapedRinging.size(); i++)
00633 m_shapedRinging[i] /= m_shapingSum;
00634
00635
00636 {
00637 std::stringstream output;
00638 output << "pmtPulse = [ ";
00639 for(int i=0; i<nPulseSamples; i++){
00640 if(i!=0) output << ", ";
00641 output << m_pmtPulse[i];
00642 }
00643 output << " ]";
00644 verbose() << output.str() << endreq;
00645 }
00646
00647 {
00648 std::stringstream output;
00649 output << "Esum Shaping Pulse = [ ";
00650 for(int i=0; i<nPulseSamples; i++){
00651 if(i!=0) output << ", ";
00652 output << m_esumResponse[i];
00653 }
00654 output << " ]";
00655 verbose() << output.str() << endreq;
00656 }
00657
00658 {
00659 std::stringstream output;
00660 output << "ringing = [ ";
00661 for(int i=0; i<nRingingSamples; i++){
00662 if(i!=0) output << ", ";
00663 output << m_ringing[i];
00664 }
00665 output << " ]";
00666 verbose() << output.str() << endreq;
00667 }
00668 {
00669 std::stringstream output;
00670 output << "shapedPmtPulse = [ ";
00671 for(int i=0; i<nPulseSamples; i++){
00672 if(i!=0) output << ", ";
00673 output << m_shapedPmtPulse[i];
00674 }
00675 output << " ]";
00676 verbose() << output.str() << endreq;
00677 }
00678 {
00679 std::stringstream output;
00680 output << "shapedRinging = [ ";
00681 for(int i=0; i<nRingingSamples; i++){
00682 if(i!=0) output << ", ";
00683 output << m_shapedRinging[i];
00684 }
00685 output << " ]";
00686 verbose() << output.str() << endreq;
00687 }
00688
00689 return StatusCode::SUCCESS;
00690 }
00691
00692 StatusCode EsIdealFeeTool::mapPulsesByChannel(
00693 const DayaBay::ElecPulseCollection::PulseContainer& pulses,
00694 std::map<DayaBay::ElecChannelId, std::vector<DayaBay::ElecPmtPulse*> >& pulseMap,
00695 int* hitCountSlices)
00696 {
00697 debug() << "Processing " << pulses.size() << " pmt pulses." << endreq;
00698 DayaBay::ElecPulseCollection::PulseContainer::const_iterator it,
00699 pulseDone = pulses.end();
00700 int timeSliceNo;
00701
00702 for (it=pulses.begin(); it != pulseDone; ++it) {
00703 DayaBay::ElecPmtPulse* pulse = dynamic_cast<DayaBay::ElecPmtPulse*>(*it);
00704 if(pulse == NULL){
00705
00706 error() << "mapPulsesByChannel(): bad pulse." << endreq;
00707 return StatusCode::FAILURE;
00708 }
00709 (pulseMap[pulse->channelId()]).push_back(pulse);
00710 timeSliceNo = int( pulse->time()/150. - 0.5 );
00711 hitCountSlices[timeSliceNo]+=int(pulse->amplitude()+0.5);
00712 }
00713 return StatusCode::SUCCESS;
00714 }
00715
00716 StatusCode EsIdealFeeTool::updateBoardSignals(DayaBay::ElecFeeCrate* crate)
00717 {
00718
00719
00720
00721 DayaBay::ElecFeeCrate::DigitalMap& nhitMap = crate->nhit();
00722 DayaBay::ElecFeeCrate::AnalogMap& esumMap = crate->esum();
00723
00724
00725 const DayaBay::ElecFeeCrate::ChannelData& channelData = crate->channelData();
00726 DayaBay::ElecFeeCrate::ChannelData::const_iterator channelIter,
00727 done = channelData.end();
00728 for(channelIter = channelData.begin(); channelIter != done; ++channelIter){
00729 DayaBay::FeeChannelId channelId = channelIter->first;
00730 const DayaBay::ElecFeeChannel& channel = channelIter->second;
00731 const DayaBay::DigitalSignal& hit = channel.hit();
00732 const DayaBay::AnalogSignal& energy = channel.energy();
00733 DayaBay::DigitalSignal& boardNhit = nhitMap[channelId.boardId()];
00734 DayaBay::AnalogSignal& boardEsum = esumMap[channelId.boardId()];
00735
00736 unsigned int hitSamples = hit.size();
00737
00738
00739
00740 if( m_hitSync.size() != hitSamples ) m_hitSync.resize( hitSamples );
00741 const int* hitStart = &hit[0];
00742 int* hitSyncStart = &m_hitSync[0];
00743 *hitSyncStart = 0;
00744 for (unsigned int i=1; i< hitSamples; i++){
00745 int* curHitSync = hitSyncStart + i;
00746 int* lastHitSync = curHitSync - 1;
00747 const int* curHit = hitStart + i;
00748 const int* lastHit = curHit - 1;
00749 if( *lastHitSync==0 && *lastHit==0 && *curHit==1 ) *curHitSync = 1;
00750 else *curHitSync = 0;
00751 }
00752
00753
00754
00755
00756
00757 if( m_hitHold.size() != hitSamples ) m_hitHold.resize( hitSamples );
00758 int* hitHoldStart = &m_hitHold[0];
00759 for (unsigned int i=0; i< hitSamples; i++) *(hitHoldStart + i) = 0;
00760
00761 if( boardNhit.size() != hitSamples ){
00762 boardNhit.resize( hitSamples );
00763 int* boardNhitStart = &boardNhit[0];
00764 for (unsigned int i=0; i< hitSamples; i++) *(boardNhitStart + i) = 0;
00765 }
00766
00767 for( unsigned int idx = 0; idx < hitSamples; idx++ ){
00768 unsigned int tempIdx = 0;
00769 if( *(hitSyncStart+idx) != 0 ){
00770 for( int holdIdx = 0; holdIdx < m_triggerWindowCycles; holdIdx++ ){
00771 unsigned int currentIdx = idx + holdIdx;
00772 if( currentIdx < hitSamples ) *(hitHoldStart + currentIdx) = 1;
00773 tempIdx++;
00774 }
00775 if(tempIdx !=0) idx = idx +tempIdx-1 ;
00776 }
00777 }
00778
00779 {
00780 int* boardNhitStart = &boardNhit[0];
00781 for(unsigned int idx = 0; idx < hitSamples; idx++)
00782 *(boardNhitStart + idx) += *(hitHoldStart + idx);
00783 }
00784 verbose() << "boardNhit " << boardNhit << endreq;
00785
00786
00787 unsigned int esumSamples = energy.size();
00788 if( boardEsum.size() != esumSamples ){
00789 boardEsum.resize( esumSamples );
00790 double* boardEsumStart = &boardEsum[0];
00791 for (unsigned int i=0; i< esumSamples; i++) *(boardEsumStart + i) = 0;
00792 }
00793 const double* energyStart = &energy[0];
00794 double* boardEsumStart = &boardEsum[0];
00795 for( unsigned int idx = 0; idx < esumSamples; idx++ ){
00796 *(boardEsumStart + idx) += *(energyStart + idx);
00797 }
00798
00799
00800
00801
00802 }
00803
00804 return StatusCode::SUCCESS;
00805 }
00806
00807 StatusCode EsIdealFeeTool::updateEsumSignals(DayaBay::ElecFeeCrate* crate){
00808
00809 DayaBay::ElecFeeCrate::AnalogMap& esumMap = crate->esum();
00810 DayaBay::ElecFeeCrate::AnalogMap::iterator boardIt = esumMap.begin();
00811
00812
00813 m_esumH.clear();
00814 m_esumL.clear();
00815 m_esumTotal.clear();
00816
00817 for (;boardIt != esumMap.end(); ++boardIt)
00818 {
00819 DayaBay::FeeChannelId boardId = boardIt->first;
00820 DayaBay::AnalogSignal& boardEnergy = boardIt->second;
00821 int boardNum = boardId.board();
00822
00823 if( m_esumL.size() != boardEnergy.size() ){
00824 m_esumL.resize( boardEnergy.size() );
00825 }
00826
00827 if( m_esumH.size() != boardEnergy.size() ){
00828 m_esumH.resize( boardEnergy.size() );
00829 }
00830
00831 if( m_esumTotal.size() != boardEnergy.size() ){
00832 m_esumTotal.resize( boardEnergy.size() );
00833 }
00834
00835 debug() << "board number " << boardNum << endreq;
00836 debug() << "Low Esum size " << m_esumL.size() << endreq;
00837 debug() << "board Energy size " << boardEnergy.size() << endreq;
00838
00839 int boardIdx = boardEnergy.size();
00840 double* boardEsumStart = &boardEnergy[0];
00841 double* esumHStart = &m_esumH[0];
00842 double* esumLStart = &m_esumL[0];
00843 double* esumTtart = &m_esumTotal[0];
00844
00845 for( int i = 0; i < boardIdx; i++)
00846 {
00847
00848 if (boardEsumType(boardNum) == DayaBay::ESumComp::kESumHigh){
00849 *(esumHStart+i) += *(boardEsumStart+i);
00850 *(esumTtart+i) += *(boardEsumStart+i);
00851
00852
00853 }
00854
00855
00856 if (boardEsumType(boardNum) == DayaBay::ESumComp::kESumLow){
00857 *(esumLStart+i) += *(boardEsumStart+i);
00858 *(esumTtart+i) += *(boardEsumStart+i);
00859
00860
00861 }
00862 }
00863
00864
00865 }
00866
00867
00868 verbose() << "Lower ESum " << m_esumL << endreq;
00869 verbose() << "Total ESum " << m_esumTotal << endreq;
00870
00871
00872 if(m_enableESumH) {
00873 verbose() << "Upper ESum " << m_esumH << endreq;
00874 shape( m_esumH, m_shapedEsumH, m_esumResponse, m_esumConvolutionWindows );
00875 m_shapedEsumH.resize( m_esumH.size() );
00876
00877
00878
00879
00880
00881 double* esumHStart = &m_shapedEsumH[0];
00882
00883 for(unsigned int i=0; i < m_shapedEsumH.size(); i++)
00884 {
00885 *(esumHStart+i) /= m_esumShapingSum;
00886
00887 }
00888
00889
00890 sample( m_shapedEsumH, crate->esumUpper(),
00891 m_simFrequency, m_eSumFrequency );
00892
00893 verbose() << "Shaped upper sum " << crate->esumUpper() << endreq;
00894 }
00895 else {
00896 m_shapedEsumH.resize( m_esumH.size() );
00897 crate->setEsumUpper(m_shapedEsumH);
00898 }
00899
00900 if(m_enableESumL) {
00901 shape( m_esumL, m_shapedEsumL, m_esumResponse, m_esumConvolutionWindows );
00902 m_shapedEsumL.resize( m_esumL.size() );
00903
00904
00905 double* esumLStart = &m_shapedEsumL[0];
00906
00907 for(unsigned int i=0; i < m_shapedEsumL.size(); i++)
00908 {
00909 *(esumLStart+i) /= m_esumShapingSum;
00910
00911 }
00912
00913
00914 sample( m_shapedEsumL, crate->esumLower(),
00915 m_simFrequency, m_eSumFrequency );
00916
00917 verbose() << "Shaped lower sum " << crate->esumLower() << endreq;
00918 }
00919 else {
00920 m_shapedEsumL.resize( m_esumL.size() );
00921 crate->setEsumLower(m_shapedEsumL);
00922 }
00923
00924 if(m_enableESumTotal) {
00925 shape( m_esumTotal, m_shapedEsumTotal, m_esumResponse, m_esumConvolutionWindows );
00926 m_shapedEsumTotal.resize( m_esumTotal.size() );
00927
00928
00929 double* esumTstart = &m_shapedEsumTotal[0];
00930
00931 for(unsigned int i=0; i < m_shapedEsumTotal.size(); i++)
00932 {
00933 *(esumTstart+i) /= m_esumShapingSum;
00934
00935 }
00936
00937
00938 sample( m_shapedEsumTotal, crate->esumTotal(),
00939 m_simFrequency, m_eSumFrequency );
00940
00941 verbose() << "Shaped total sum " << crate->esumTotal() << endreq;
00942 }
00943 else {
00944 m_shapedEsumTotal.resize( m_esumTotal.size() );
00945 crate->setEsumTotal(m_shapedEsumTotal);
00946 }
00947
00948
00949
00950 digitize(crate->esumTotal(), crate->esumADC(), m_adcRange, m_adcBits, m_adcOffset);
00951
00952 verbose() << "Digitized signal " << crate->esumADC() << endreq;
00953
00954 return StatusCode::SUCCESS;
00955 }
00956
00957
00958 DayaBay::ESumComp::ESumComp_t EsIdealFeeTool::boardEsumType(int boardId){
00959
00960
00961 if (boardId >= 1 && boardId < 7)
00962 {
00963 return DayaBay::ESumComp::kESumLow;
00964 }
00965
00966 if (boardId >= 7 && boardId < 13)
00967 {
00968 return DayaBay::ESumComp::kESumHigh;
00969 }
00970
00971 if (boardId == 13)
00972 {
00973 return DayaBay::ESumComp::kESumNone;
00974 }
00975
00976 return DayaBay::ESumComp::kUnknown;
00977 }
00978
00979
00980 void EsIdealFeeTool::convolve(const DayaBay::AnalogSignal& signal,
00981 const DayaBay::AnalogSignal& response,
00982 DayaBay::AnalogSignal& output)
00983 {
00984
00985
00986
00987 int sN = signal.size();
00988 int rN = response.size();
00989 output.resize(sN+rN);
00990 const double* sigD = &(signal[0]);
00991 const double* resD = &(response[0]);
00992 double* outD = &(output[0]);
00993 for (int i=0; i<(sN+rN); i++) {
00994 double sum = 0;
00995 for (int j=0; j<rN; j++) {
00996 if (resD[j]==0)
00997 continue;
00998 if (i-j>=0 && i-j<sN)
00999 sum+=sigD[i-j]*resD[j];
01000 }
01001 outD[i]=sum;
01002 }
01003 }
01004
01005
01006 void EsIdealFeeTool::shape(const DayaBay::AnalogSignal& rawSignal,
01007 DayaBay::AnalogSignal& shapedSignal,
01008 const DayaBay::AnalogSignal& response, std::map<int,int> convolutionWindows){
01009 std::map<int,int>::iterator convIt = convolutionWindows.begin();
01010 int convStart, convEnd;
01011 shapedSignal.clear();
01012 int samples = rawSignal.size();
01013 shapedSignal.resize(samples, 0);
01014 DayaBay::AnalogSignal rawFrag;
01015 DayaBay::AnalogSignal shapedFrag;
01016 for(; convIt != convolutionWindows.end(); convIt++){
01017 convStart = convIt->first * 96;
01018 convEnd = convIt->second * 96;
01019 if (convEnd > samples) convEnd = samples;
01020 rawFrag.resize( (convEnd - convStart) );
01021 shapedFrag.resize( (convEnd - convStart) );
01022 std::copy( rawSignal.begin() + convStart, rawSignal.begin() + convEnd , rawFrag.begin() );
01023 convolve( rawFrag, response, shapedFrag );
01024 shapedFrag.resize( (convEnd - convStart) );
01025 std::copy( shapedFrag.begin(), shapedFrag.end(), shapedSignal.begin() + convStart );
01026 debug() << "Convolve from " << convStart *1.5625 << " ns to " << convEnd * 1.5625 << " ns" << endreq;
01027 }
01028 }
01029 void EsIdealFeeTool::sample(const DayaBay::AnalogSignal& signal,
01030 DayaBay::AnalogSignal& output,
01031 int inputFrequency,
01032 int outputFrequency)
01033 {
01034
01035
01036 if(inputFrequency == outputFrequency){
01037 if( output.size() != signal.size() ) output.resize( signal.size() );
01038 memcpy(&output[0],&signal[0],output.size()*sizeof(double));
01039 return;
01040 }
01041 unsigned int nSamples = convertClock(signal.size(),
01042 inputFrequency,
01043 outputFrequency);
01044 if(output.size() != nSamples){
01045 if(output.capacity() > nSamples){
01046 output.resize(nSamples);
01047 }else{
01048
01049 DayaBay::AnalogSignal* newSignal = new DayaBay::AnalogSignal(nSamples);
01050 output.swap(*newSignal);
01051 delete newSignal;
01052 }
01053 }
01054
01055 if(inputFrequency < outputFrequency){
01056 warning() << "sample(): Input frequency " << inputFrequency
01057 << " is less than the sampling frequency " << outputFrequency
01058 << endreq;
01059 }
01060 if( (inputFrequency % outputFrequency) == 0 ){
01061
01062 int relativeFrequency = inputFrequency / outputFrequency;
01063 double* outputStart = &output[0];
01064 const double* inputStart = &signal[0];
01065 int inputIdx = 0;
01066 unsigned int outputIdx = 0;
01067 while(outputIdx != nSamples){
01068 *(outputStart + outputIdx) = *(inputStart + inputIdx);
01069 outputIdx++;
01070 inputIdx+=relativeFrequency;
01071 }
01072 return;
01073 }else{
01074
01075 double* outputStart = &output[0];
01076 const double* inputStart = &signal[0];
01077 int inputIdx = 0;
01078 unsigned int outputIdx = 0;
01079 double freqScale = inputFrequency / outputFrequency;
01080 while( outputIdx != nSamples ){
01081 inputIdx = int(outputIdx * freqScale);
01082 *(outputStart + outputIdx) = *(inputStart + inputIdx);
01083 outputIdx++;
01084 }
01085 }
01086 }
01087
01088 void EsIdealFeeTool::digitize(const DayaBay::AnalogSignal& signal,
01089 DayaBay::DigitalSignal& output,
01090 double range, int bits, double offset)
01091 {
01092
01093 if(output.size() != signal.size()) output.resize(signal.size());
01094 int maxADC = int( pow(2, bits) );
01095 int signalN = signal.size();
01096 const double* inputStart = &signal[0];
01097 int* outputStart = &output[0];
01098 for (int i=0; i<signalN; i++) {
01099 int* curOutput = outputStart + i;
01100 *(curOutput) = int(((*(inputStart + i))/range)*maxADC + offset);
01101 if (*curOutput>maxADC)
01102 *curOutput=maxADC;
01103 if (*curOutput<0)
01104 *curOutput=0;
01105
01106
01107
01108
01109
01110
01111
01112 }
01113 }
01114
01115 void EsIdealFeeTool::discriminate(const DayaBay::AnalogSignal& signal,
01116 double threshold,
01117 std::vector<int>& output,
01118 DayaBay::Threshold::Threshold_t mode)
01119 {
01120
01121
01122
01123
01124
01125 threshold *= m_speAmp;
01126 threshold *= m_gainFactor;
01127 threshold *= m_discThreshScale;
01128
01129 output.clear();
01130 const double* signalStart = &signal[0];
01131 unsigned int nSamples = signal.size();
01132 if (DayaBay::Threshold::kCrossing == mode) {
01133 bool aboveThreshold = false;
01134 for (unsigned int i=0;
01135 i < nSamples;
01136 ++i) {
01137 if (!aboveThreshold && (*(signalStart+i))>=threshold) {
01138 output.push_back(i);
01139 aboveThreshold = true;
01140 } else if (aboveThreshold && (*(signalStart+i))<threshold)
01141 aboveThreshold=false;
01142 }
01143 } else if (DayaBay::Threshold::kAbove == mode) {
01144 for (unsigned int i=0;
01145 i < nSamples;
01146 ++i)
01147 if ((*(signalStart+i))>=threshold)
01148 output.push_back(i);
01149 } else {
01150 error() << "discriminate(): Unknown Discrimination Mode" << endreq;
01151 }
01152 }
01153
01154 double EsIdealFeeTool::pmtPulse(double deltaT, int nPulse, bool addOvershoot) {
01155
01156 double shift = -2.5e-9;
01157 if (deltaT-shift<0) return 0.;
01158 double v0 = -m_speAmp;
01159 double width;
01160 double mu;
01161 if ( nPulse != 0 ){
01162 width = getPulseWidth( nPulse);
01163 mu = getPulseForm( nPulse);
01164 }
01165 else {
01166 width = 7.5e-9;
01167 mu = 0.5;
01168 }
01169 if(deltaT > 100 * width) return 0;
01170 double amplitude = 0;
01171 if(deltaT < 10 * width){
01172
01173 amplitude += m_gainFactor* v0 * exp( -pow( log( (deltaT-shift)/width),2)
01174 / ( 2*pow(mu,2) ) ) * Gaudi::Units::volt;
01175 }
01176 if( addOvershoot ){
01177
01178 amplitude += overshoot( deltaT, nPulse );
01179 }
01180 return amplitude;
01181 }
01182
01183 double EsIdealFeeTool::pulseShaping(double deltaT) {
01184
01185
01186
01187
01188 if (deltaT<0)
01189 return 0.;
01190 double t0 = 25.0e-9;
01191 double r = deltaT/t0;
01192 return pow( r, 4 ) * exp( -r );
01193 }
01194
01195 double EsIdealFeeTool::esumShaping(double deltaT) {
01196
01197 if (deltaT<0)
01198 return 0.;
01199 double t0 = 56.0;
01200 double r = deltaT/t0;
01201 return exp( -r )*(1/t0);
01202 }
01203
01204 double EsIdealFeeTool::noise() {
01205
01206
01207
01208 return m_gauss() * m_noiseAmp;
01209 }
01210
01211 int EsIdealFeeTool::convertClock(int clock, int inputFrequency, int outputFrequency){
01212
01213
01214 return int( ( clock / double(inputFrequency) ) * outputFrequency );
01215 }
01216
01217 double EsIdealFeeTool::satFactor(double amp){
01218
01219 if (amp<0) return 1.;
01220 double y = 1.e6 * amp;
01221 if(y < 1)
01222 return 1. - 0.5311049 * y + 0.1818968 * pow(y,2);
01223 if(y < 9.5)
01224 return 0.778264 - 0.1400074 * y + 0.01303484 * pow(y,2) - 0.0004701175 * pow(y,3) + 0.000002450386 * pow(y,4);
01225 return 0.4212711 - 2.560691e-02 * y + 8.003103e-04 * pow(y,2) - 1.226207e-05 * pow(y,3) + 7.282761e-08 * pow(y,4);
01226 }
01227
01228 double EsIdealFeeTool::overshoot(double deltaT, int nPulse) {
01229
01230 if (deltaT < 0) return 0.;
01231 double v0 = getOvershootAmp( nPulse );
01232 double tau = 155.e-9;
01233 double expo = v0*exp(-deltaT/tau);
01234 double t0 = 50e-9;
01235 double t1 = 10e-9;
01236 double fermi = 1. / (exp( (t0 - deltaT) / t1) + 1.);
01237 return fermi * expo * Gaudi::Units::volt;
01238 }
01239
01240 double EsIdealFeeTool::ringing(double deltaT, double pulseMax) {
01241
01242 if (deltaT<0) return 0.;
01243 double v0 = 0.00154 * pulseMax;
01244 double period = 1.63e-6;
01245 double shift= 0.99e-6;
01246 double expo = exp(-deltaT/12.e-6);
01247 double t0 = 0.6e-6;
01248 double t1 = 0.1e-6;
01249 double fermi = 1./(exp((t0-deltaT)/t1)+1);
01250 return v0 * fermi * expo * cos(6.28*(deltaT-shift)/period);
01251 }
01252
01253 double EsIdealFeeTool::getPulseWidth(int nPulse){
01254
01255 if( nPulse > 2000) return 20.7842e-9;
01256 double x = double(nPulse);
01257 double t0 = 1200.;
01258 double t1 = 200.;
01259 double f = 1. / (exp( (t0-x)/t1 ) + 1.);
01260 return (7.+12.*f+x/1000.)*1.e-9;
01261 }
01262
01263 double EsIdealFeeTool::getPulseForm(int nPulse){
01264
01265 if( nPulse > 3000) return 0.2;
01266 double x = double(nPulse);
01267 double t0 = 900.;
01268 double t1 = 200.;
01269 double f = 0.45 - 0.25 * 1. / (exp( (t0-x)/t1) + 1.);
01270 return f;
01271 }
01272
01273 double EsIdealFeeTool::getOvershootAmp(int nPulse){
01274
01275 if( nPulse > 2000) return 2.4033e-4;
01276 double x = double(nPulse);
01277 double t0 = -0.0148227;
01278 double t1 = -1186.36;
01279 double t2 = 287.821;
01280 double t3 = 0.0150628;
01281 double f = t0 / (exp( (t1-x)/t2) + 1.) + t3;
01282 return f;
01283 }