SimVaElectronics.cxx

Go to the documentation of this file.

#include <cmath>
#include <TMath.h>
#include "MessageService/MsgService.h"
#include "Digitization/DigiSignal.h"
#include "SimVaElectronics.h"
#include "SimElecMaker.h"
#include "Calibrator/Calibrator.h"
#include "Calibrator/CalVaLinearity.h"

// Register this class with the factory.
SimElecMakerProxy<SimVaElectronics> gSimVaElec("SimVaElectronics");

CVSID("$Id: SimVaElectronics.cxx,v 1.30 2008/10/24 12:01:43 blake Exp $");

ClassImp(SimVaElectronics)

SimVaElectronics::SimVaElectronics(  VldContext context, TRandom* random )
  : SimElectronics( context, random ),
    fVaGain(0.375/Munits::fC),
    fVaPedWidth(2.8),
    fTimeSmearingMode(SimTimeSmearingMode::kBreitWigner),
    fTimingWidth(0.0),
    fVaSparsifyThresh(17),
    fVaDynodeThresh(57),
    fVaNonlinearity(1),
    fVaLinearityP0(24),
    fVaLinearityP1(16091.3),
    fVaLinearityP2(4.63),
    fVaLinearityP3(13683.9),
    fVaLinearityP4(6650.11),
    fVarcTriggerMode(SimVarcTriggerMode::k2of36),
    fVarcTriggerWindow(400.*Munits::ns),
    fVaChipRandomDeadRate(700*Munits::hertz),
    fVaChipDeadTime(5.*Munits::microsecond),
    fDoLookupNonlinearity(0),
    fCalVaLinearity(0)
{
  fCalVaLinearity = new CalVaLinearity(RawChannelId(),
                                       fVaLinearityP0, fVaLinearityP1,
                                       fVaLinearityP2, fVaLinearityP3,
                                       fVaLinearityP4
                                       );
}


SimVaElectronics::~SimVaElectronics()
{
  if(fCalVaLinearity) delete fCalVaLinearity;
}


void SimVaElectronics::ReadoutDetector( SimPmtList& pmtList )
{
  PlexHandle plex(fContext);

  // For timing-cut version of varc pretrigger.
  std::map<int,std::vector<double> > varcTrigTimes;
  std::map<int,std::vector<double> >::iterator it;
  
  // For Dead Chips.
  std::map<int,bool> deadChips;
  
  // Do the dead chips:
  if(fVaChipRandomDeadRate>0.) { // don't bother if turned off
    SimPmtList::iterator PmtIt;
    for(PmtIt = pmtList.begin(); PmtIt!= pmtList.end(); PmtIt++)  {
      // how likely is this pmt to be dead at this time?
      // Ignores finite event length
      double deadProb = fVaChipRandomDeadRate*fVaChipDeadTime;
      // Roll the dice:
      bool isDead = false;
      if(fRandom->Rndm()<deadProb) isDead = true;
      // Remember this result:
      deadChips[PmtIt->first] = isDead;
    }
  }

  Int_t  crate, varc, vmm;
  
  if(fVarcTriggerMode != SimVarcTriggerMode::kNone) {    
    // First pass: find dynode trigger times
    
    SimPmtList::iterator PmtIt;
    for(PmtIt = pmtList.begin(); PmtIt!= pmtList.end(); PmtIt++)  {
      SimPmt* pmt = PmtIt->second;
      if(pmt) {
        // Make sure the PMT reads out to VA electronics:
        RawChannelId rcid = plex.GetRawChannelId(pmt->GetPixelSpotId(1));
        if ( (rcid.GetElecType()==ElecType::kVA) ) {

          // Apply dynode trigger; don't read out if below thresh.
          // Don't read out if dead chip, either.
          if((DynodeTrigger(pmt->GetDynodeCharge())) && (!deadChips[PmtIt->first])) {
            crate = rcid.GetCrate();
            varc  = rcid.GetVarcId();
            vmm   = rcid.GetVmm();
            
            // Compose a unique index for this varc (or vmm)
            int index = crate*100 + varc*10;
            if(fVarcTriggerMode == SimVarcTriggerMode::k2of6) {
              index += vmm;
            }
            
            Double_t time = pmt->GetDynodeTime();
            varcTrigTimes[index].push_back(time);
          }
        }      
      }
    }
        
    // For each varc or vmm:
    for(it = varcTrigTimes.begin(); it!=varcTrigTimes.end(); it++) {
      sort(it->second.begin(),it->second.end()); // Sort the trigger times.
      // Debugging output...
      //for(UInt_t i=0;i<it->second.size();i++) {
      //        cout << it->first << "\t" 
      //             << it->second[i]*1e9 << endl;
      //}
    }    
  }
  

  // Apply pretrigger and read out
  SimPmtList::iterator PmtIt;
  for(PmtIt = pmtList.begin(); PmtIt!= pmtList.end(); PmtIt++)  {
    SimPmt* pmt = PmtIt->second;
    if(pmt) {
      // Make sure the PMT reads out to VA electronics:
      RawChannelId rcid = plex.GetRawChannelId(pmt->GetPixelSpotId(1));
      if ( rcid.GetElecType()==ElecType::kVA ) {
        // Apply dynode trigger; don't read out if below thresh.
        crate = rcid.GetCrate();
        varc  = rcid.GetVarcId();
        vmm   = rcid.GetVmm();
        
        bool trig = true;

        if(fVarcTriggerMode != SimVarcTriggerMode::kNone) {
          // Compose a unique index for this varc (or vmm)
          int index = crate*100 + varc*10;
          if(fVarcTriggerMode == SimVarcTriggerMode::k2of6) {
            index += vmm;
          }
          
          // Assume it's bad
          trig = false;
          
          std::vector<double> &tVec = varcTrigTimes[index];
          // Look for another hit within the trigger window.
          int nhits = 0;
          for(UInt_t i=0; i<tVec.size(); i++) {
            double dt = tVec[i] - pmt->GetDynodeTime();
            if(fabs(dt)<fVarcTriggerWindow)  // There's another trigger in the window
              nhits++;
          }
          if(nhits>1) trig=true; // Don't count itself
        }
        
        if(DynodeTrigger(pmt->GetDynodeCharge()))
          ReadoutPmt(pmt,trig, deadChips[PmtIt->first]);          
      }
    }
  }
  
}



void SimVaElectronics::ReadoutPmt( SimPmt* pmt, bool trig, bool isDead)
{
  // Reads out a single PMT.
  // Adds all readout to the list of SimDigits.

  MSG("DetSim",Msg::kDebug) << "SimVaElectronics::ReadoutPmt " 
                            << pmt->GetTubeId().AsString() 
                            << " type " << pmt->GetType()
                            << endl;

  PlexHandle plex(fContext);

  // Interesting statistics:
  if((trig)&&(!isDead)) {
    AddDigitsAfterFETrigger( pmt->GetTotalHitPixels(true) );
    AddAdcsAfterFETrigger( pmt->GetTotalCharge() * fVaGain );
  }

  for(int ipixel = 1; ipixel<= pmt->GetNumberOfPixels(); ipixel++) {
    RawChannelId rcid = plex.GetRawChannelId(pmt->GetPixelSpotId(ipixel));

    // Find which spot had the biggest charge.
    // We use this spot to look up the nonlinearity calibration.
    Float_t bigpe=0;
    Int_t   bigspot=1;
    for(int ispot=1;ispot<=pmt->GetNumberOfSpots();ispot++) {
      float pe = pmt->GetPeXtalk(ipixel, ispot);
      if(pe > bigpe) {
        bigpe = pe;
        bigspot = ispot;
      }
    }

    DigiSignal* signal = pmt->CreateSignal(ipixel);
    AddSignal(signal);
    int errorbits = 0;
    if(!trig) errorbits+=1;   // Pretrigger failed
    if(isDead) errorbits+=2;  // Dead chip

    if(!(rcid.IsNull())) {
      Int_t adc = GenSimulatedADC(pmt->GetCharge(ipixel),
                                  rcid, 
                                  pmt->GetPixelSpotId(ipixel, bigspot) );
      Int_t tdc = GenSimulatedTDC(pmt->GetDynodeTime(), rcid);
      
      SimDigit d(  pmt->GetPixelSpotId(ipixel),
                   rcid,
                   adc,
                   tdc,
                   signal,
                   errorbits
                   );
      MSG("DetSim",Msg::kVerbose) << "ReadoutPmt: " << d.AsString() << std::endl;    
      
      // Sparsify and record.
      if(d.GetADC() > fVaSparsifyThresh) { 
        AddDigit(d);
        
        // Stats:
        AddDigitsAfterSpars(1);
        AddAdcsAfterSpars(d.GetADC());
      }
    }
  }
}


Bool_t SimVaElectronics::DynodeTrigger( double dynCharge )
{
  // Returns true if the PMT would fire the dynode trigger.
  
  // ASDLite DAC setttings: (snip email by Roy Lee)
  //The DAC setting for FD dynode thresholds go from 63 (no threshold) to 0
  //(maximum thresold), where the maximum corresponds to 160 fC +- 10%,
  //assuming linearity over full scale.

  // According to Alfons, from dynode trigger scans, 1 dynode trigger count ~= 2.4 ADC counts.
  // Remember that the charge on the dynode is 0.8 of the charge on the anode, so 
  // the the maximum corresponds to 63 * 2.4 * 0.8 ADCs.

  // These two numbers seem to disagree. This is resolved because Roy forgot to mention that the
  // dynode charge is split into TWO capacitors on the base, so the trigger level of the ASDLite
  // is effectively 1/2 of what he said.

  // Thus, the full-scale trigger threshold is 322 fC (dynode) or ~160 fC (on the asdlite cap)
  
  // Threshold for max settting:
  const float kThresh0 = 63.0 * 2.4 *0.8 /(fVaGain); // fC
  
  double thresh = (double)(63 - fVaDynodeThresh)/63. * kThresh0;

  if(dynCharge > thresh) return true;
  else return false;
}

int SimVaElectronics::GenSimulatedADC( double inCharge, 
                                       const RawChannelId&,
                                       const PlexPixelSpotId& psid)
{
  // Apply nonlinearity.
  double adc = inCharge * fVaGain;

  if(fDoLookupNonlinearity) {
    PlexHandle plex(fContext);
    PlexStripEndId seid = plex.GetStripEndId(psid);
    if(seid.IsValid()) {
      adc = Calibrator::Instance().DecalLinearity(adc, seid);
    }
  }
  else if(fVaNonlinearity) {
    if(adc > 16384) adc = 16384; // Truncate to region where UnLinearize doesn't work.
    adc = fCalVaLinearity->UnLinearize(adc);
  }

  double x = fRandom->Gaus(adc, fVaPedWidth);
  return (int)(TMath::Floor(x));
}

int SimVaElectronics::GenSimulatedTDC( double inTime, const RawChannelId& rcid  )
{
  
  // VA electronics, tick = 1.5625 ns (640MHz)
  if(inTime>=kPmtTime_Never) {
    MSG("DetSim",Msg::kWarning) << "GenSimulatedTDC got a hit with a PMT that never fired!" << endl;
    // Return a false number so this hit doesn't correlate with anything real.
    return 650000000;
  }

  // smear the raw time by fTimingWidth (in nanoseconds)
  double deltaTime = 0.0;
  if( fTimeSmearingMode>0 && fTimingWidth>0.0 ) {

    // smearing with Breit-Wigner function
    if( fTimeSmearingMode==SimTimeSmearingMode::kBreitWigner ){
      deltaTime = ( fRandom->BreitWigner( 0.0, fTimingWidth ) )*Munits::ns;
    }

    // smearing with Gaussian function
    else if( fTimeSmearingMode==SimTimeSmearingMode::kGaussian ){
      deltaTime = ( fRandom->Gaus( 0.0, fTimingWidth ) )*Munits::ns;
    }

    // smearing with Uniform function
    else if( fTimeSmearingMode==SimTimeSmearingMode::kUniform ){
      deltaTime = ( fRandom->Uniform( -fTimingWidth, +fTimingWidth ) )*Munits::ns;
    }
    
  }

  //const double tdc_convert = 1.0/1.5625e-9;
  double x = Calibrator::Instance().GetTDCFromTime(inTime+deltaTime,rcid);

  return (int)(TMath::Floor(x));
}

void
SimVaElectronics::Print(Option_t*) const
{
  printf("VA Front End:  vaGain               %f ADC/fC\n",fVaGain*Munits::fC);
  printf("               vaPedWidth           %f ADCs\n",fVaPedWidth);
  printf("               vaTimingWidth        %f ns\n",fTimingWidth);
  printf("               vaSparsifyThresh     %d ADCs\n",fVaSparsifyThresh);
  printf("               vaDynodeThresh       %d DACs\n",fVaDynodeThresh);
  printf("               vaNoninearity        %d (%s)\n",fVaNonlinearity,(fVaNonlinearity)?"Nonlinear":"Linear");
  printf("               vaLinearityParams:   %f %.1f %.1f %.1f %.1f\n",fVaLinearityP0, fVaLinearityP1,fVaLinearityP2,fVaLinearityP3,fVaLinearityP4);
  printf("VA Back End:   varcTriggerMode      %d (0=none, 1=2/6, 2=2/36)\n",fVarcTriggerMode);
  printf("               varcTriggerWindow    %f ns\n",fVarcTriggerWindow);
  printf("               vaChipRandomDeadRate %f Hz\n",fVaChipRandomDeadRate);
  printf("               vaChipDeadTime       %f us\n",fVaChipDeadTime/Munits::microsecond);
  printf("               doLookupNonlinearity  %s \n",fDoLookupNonlinearity?("on"):("off") );

} 
 
void   
SimVaElectronics::Config( Registry& config )
{
  // Modify the configuration.
  //
  // Use this method to set static members with the class configuration.
  config.Get("vaGain",fVaGain);
  config.Get("vaPedWidth",fVaPedWidth);
  config.Get("vaTimeSmearingMode",fTimeSmearingMode);
  config.Get("vaTimingWidth",fTimingWidth);
  config.Get("vaSparsifyThresh",fVaSparsifyThresh);
  config.Get("vaDynodeThresh",fVaDynodeThresh);
  config.Get("vaNonlinearity",fVaNonlinearity);
  config.Get("vaLinearityP0",fVaLinearityP0);
  config.Get("vaLinearityP1",fVaLinearityP1);
  config.Get("vaLinearityP2",fVaLinearityP2);
  config.Get("vaLinearityP3",fVaLinearityP3);
  config.Get("vaLinearityP4",fVaLinearityP4);  
  config.Get("varcTriggerMode",fVarcTriggerMode);
  config.Get("varcTriggerWindow",fVarcTriggerWindow);
  config.Get("vaChipRandomDeadRate",fVaChipRandomDeadRate);
  config.Get("vaChipDeadTime",fVaChipDeadTime);
  config.Get("doLookupNonlinearity",fDoLookupNonlinearity);
  
  if(fCalVaLinearity) delete fCalVaLinearity;
  fCalVaLinearity = new CalVaLinearity(RawChannelId(),
                                       fVaLinearityP0,
                                       fVaLinearityP1,
                                       fVaLinearityP2,
                                       fVaLinearityP3,
                                       fVaLinearityP4
                                       );

  if(fVaSparsifyThresh<=0) 
    MSG("DetSim",Msg::kWarning) << "WARNING: The VA sparsification threshold " << std::endl
                                << "         is less than zero!" << std::endl;                 
  SimElectronics::Config(config);
}

---