SimQieElectronics.cxx

Go to the documentation of this file.

#include <TMath.h>
#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "Digitization/DigiSignal.h"
#include "SimQieElectronics.h"
#include "SimElecMaker.h"
#include "Calibrator/Calibrator.h"

// Register this class with the factory.
SimElecMakerProxy<SimQieElectronics> gSimQieElec("SimQieElectronics");

CVSID("$Id: SimQieElectronics.cxx,v 1.21 2008/10/15 16:46:34 tagg Exp $");

ClassImp(SimQieElectronics)

SimQieElectronics::SimQieElectronics(  VldContext context, TRandom* random )
  : SimElectronics( context, random ),
    fClock(context),
    fQieDoLookup(1),
    fQieDacCharge(1.4*Munits::fC),
    fQieDacPerRangeZeroAdc(1.92),
    fQiePedestalWidthAdc(1.5),
    fQieAdcRms(0),
    fQiePedestalDac(50),
    fQieSparsifyThresh(70),
    fDoLookupNonlinearity(0)
{
  SetupTables();
}

void SimQieElectronics::Reset( const VldContext& newContext )
{
  fClock.Reset(newContext);
  SimElectronics::Reset(newContext);
}


void SimQieElectronics::ReadoutDetector( SimPmtList& pmtList )
{
  PlexHandle plex(fContext);
  SimPmtList::iterator it;
  for(it = pmtList.begin(); it!= pmtList.end(); it++)  {
    SimPmt* pmt = it->second;
    if(pmt) {
      RawChannelId rcid = plex.GetRawChannelId(pmt->GetPixelSpotId(1));
      // Make sure the PMT reads out to VA electronics:      
      if ( (rcid.GetElecType()==ElecType::kQIE) ) {       
        ReadoutPmt( pmt );
      }
    }
  }
  
}

void SimQieElectronics::ReadoutPmt( SimPmt* pmt )
{
  // Reads out a single PMT.
  // Adds all readout to the list of SimDigits.
  
  // This vaugely simulates the real ND,
  // which reads out all buckets from a channel before
  // moving to the next one. This is only incorrect in that
  // this uses pixel order instead of master/minder address order. 
  // Shouldn't matter.

  // Also, this assumes a fast extraction mode, such that no 
  // hits are out of the readout window.  For resonant extraction
  // or cosmics, you need to look for a dynode trigger, find latch the
  // dynode trigger time, and put a validity window around the dynode
  // trigger QIE tick.  It's irrelevant for neutrino physics (at the moment)
  // so I'm ignoring it.  --N

  // This does NOT correctly deal with multiplexing!  Must fix!

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

  PlexHandle plex(fContext);

  // No trigger, so add up all hits
  AddDigitsAfterFETrigger( pmt->GetTotalHitPixels(true) );
  AddAdcsAfterFETrigger( pmt->GetTotalCharge() * fQieDacCharge );
    
  // Build a list of channels we can read out.
  std::map<UInt_t,std::vector<UInt_t> > channels;
  std::map<UInt_t,std::vector<UInt_t> >::iterator chItr;
  
  for(int pixel = 1; pixel<= pmt->GetNumberOfPixels(); pixel++) {
    PlexPixelSpotId psid = pmt->GetPixelSpotId(pixel);
    RawChannelId    rcid = plex.GetRawChannelId(psid);
    channels[rcid.GetEncoded()].push_back(pixel);
  }
  
  // Loop through buckets.
  SimPmtBucketIterator bIt( *pmt );
  for( ; !bIt.End(); bIt.Next() ) {
    
    Int_t bucketID = bIt.BucketId();
    SimPmtTimeBucket& bucket = bIt.Bucket();

    for( chItr = channels.begin(); chItr != channels.end(); chItr++) {
      Float_t totalCharge = 0;
      DigiSignal* totalSignal = new DigiSignal();
      AddSignal(totalSignal);
      RawChannelId    rcid(chItr->first);
      PlexPixelSpotId psid;

      for(UInt_t ipix = 0; ipix < chItr->second.size(); ipix++) {
        Int_t pixel = chItr->second[ipix];
        psid = pmt->GetPixelSpotId(pixel);
        SimPixelTimeBucket& pixBucket = bucket.GetPixelBucket(pixel);

        DigiSignal* signal = pixBucket.CreateSignal();
        Float_t     charge = pixBucket.GetCharge();

        totalSignal->Merge(*signal);
        totalCharge+=charge;

        delete signal;
      }

      // Protect against invalid channels on tubes that have some valid pixels.
      if(!(rcid.IsNull())) {            
        if((totalCharge<-1000*Munits::fC)||(totalCharge>1e6*Munits::fC)) {
          MSG("DetSim",Msg::kError) << "Abnormal charge " << totalCharge/Munits::fC 
                                    << " fC in pmt " << pmt->GetTubeId().AsString() << endl;
          pmt->Print();
        } 
        SimDigit d(  psid,                        // Pixel, or one of them
                     rcid,                        // raw channel
                     GenSimulatedADC( rcid, psid, totalCharge, bucketID ), // ADC
                     bucketID,                        // The TDC _is_ the bucket ID.
                     totalSignal,                     // Truth
                     0                                // Error bits.
                     );
        MSG("DetSim",Msg::kVerbose) << "SimQieElectronics::ReadoutPmt: " << d.AsString() << std::endl;
        
        // Sparsify and record.
        if(d.GetADC() > fQieSparsifyThresh) {
          AddDigit(d);
          
          // Stats:
          AddDigitsAfterSpars(1);
          AddAdcsAfterSpars(d.GetADC());
        }
      }
    }
  }
}


Bool_t SimQieElectronics::DynodeTrigger( SimPmt*  )
{
  // No dynode simulation...yet.
  return true;  
}



void SimQieElectronics::SetupTables( void )
{
  // These are the slopes (ADC/DAC) of each range.
  fAdcDacSlope[0] = 1./fQieDacPerRangeZeroAdc; // By definition.
  // Each higher range has double the gain of the last one.
  // These values could be improved still; I believe they are just approximate.
  fAdcDacSlope[1] = fAdcDacSlope[0]*0.5;
  fAdcDacSlope[2] = fAdcDacSlope[1]*0.5;
  fAdcDacSlope[3] = fAdcDacSlope[2]*0.5;
  fAdcDacSlope[4] = fAdcDacSlope[3]*0.5;
  fAdcDacSlope[5] = fAdcDacSlope[4]*0.5;
  fAdcDacSlope[6] = fAdcDacSlope[5]*0.5;
  fAdcDacSlope[7] = fAdcDacSlope[6]*0.5;

  // These values are taken from a plot from Charlie Nelson. See
  // MENUswitchpoints.pdf in the doc directory.
  // These need not be terribly accurate, since they very from channel-to-channel
  // by 4 or 6 counts RMS.
  fFlipLowAdc[0] = 0;     // Min
  fFlipLowAdc[1] = 23.7;
  fFlipLowAdc[2] = 25.6;
  fFlipLowAdc[3] = 26.5;
  fFlipLowAdc[4] = 27.1;
  fFlipLowAdc[5] = 27.0;
  fFlipLowAdc[6] = 27.2;
  fFlipLowAdc[7] = 29;

  fFlipHighAdc[0] = 185;
  fFlipHighAdc[1] = 190;
  fFlipHighAdc[2] = 196;
  fFlipHighAdc[3] = 199;
  fFlipHighAdc[4] = 192;
  fFlipHighAdc[5] = 200;
  fFlipHighAdc[6] = 177;
  fFlipHighAdc[7] = 255; // Max

  // This is the point, in DAC counts, where the range flips.
  // i.e. fFlipPoint[1] is the charge where the QIE flips from range 0 to 1
  // This is determined entirely from the above tables.
  fFlipPoint[0] = 0;
  for(int i=1;i<9;i++)
    fFlipPoint[i] = fFlipPoint[i-1] + (fFlipHighAdc[i-1]-fFlipLowAdc[i-1])/fAdcDacSlope[i-1];
}




int SimQieElectronics::GenSimulatedADC( RawChannelId, 
                                        PlexPixelSpotId psid,
                                        double inCharge,
                                        int tdc)
{
  //
  // This routine converts from a charge in fC to a 
  // QIE digitized charge.  Note that 'ADC' is an overloaded term here,
  // since these aren't ADC counts but rather lookup table values.
  // The units are actually DAC counts, which are 1.4 fC/ DAC.
  //

  // Convert to DAC charge units. These are the units of charge
  // used by the sharge injection system.
  double true_dacs = inCharge / fQieDacCharge;

  // Do nonlinearity curve.
  if (fDoLookupNonlinearity) {
    PlexHandle plex(fContext);
    PlexStripEndId seid = plex.GetStripEndId(psid);
    if(seid.IsValid())
      true_dacs = Calibrator::Instance().DecalLinearity(true_dacs, seid);
  }

  // Add some smear, by the pedestal RMS in DAC counts.
  double dacs = fRandom->Gaus(true_dacs, 
                              fQiePedestalWidthAdc * fQieDacPerRangeZeroAdc);

  // Figure out what range you're in.
 // Find the range and pos within range.
  int range = 0;
  if     (dacs < fFlipPoint[1])  range = 0;
  else if(dacs < fFlipPoint[2])  range = 1;
  else if(dacs < fFlipPoint[3])  range = 2;
  else if(dacs < fFlipPoint[4])  range = 3;
  else if(dacs < fFlipPoint[5])  range = 4;
  else if(dacs < fFlipPoint[6])  range = 5;
  else if(dacs < fFlipPoint[7])  range = 6;
  else                           range = 7;


  // Now find the ADC value.
  double adc = fFlipLowAdc[range] + fAdcDacSlope[range]*(dacs - fFlipPoint[range]);

  // Smear by a constant resolution function.
  // i.e. if all values are wrong by 1.2 ADC counts rms, this does it.
  // This is set to zero by default because it is difficult to seperate
  // the charge injector error from the ADC conversion error.
  if(fQieAdcRms>0)
    adc = fRandom->Gaus(adc,fQieAdcRms);

  // Make it an integer again.
  // Rounding is probably the most appropriate, given that we have correctly 
  // calibrated the pedestals and whatnot.
  
  if(adc>255) adc=255.;
  if(adc<0)   adc=0.;
  double iadc = TMath::Nint(adc);

  if(fQieDoLookup==0) {
    // We have turned off the lookup table, so we need to simply pack the output
    // word:
    int adc8bits = (int)iadc;
    int capid = tdc%4; // CapID rotates 0123 with the clock
    int dataword = (adc8bits & 0xFF) | (range << 8)  | (capid << 11); 
    return dataword;
  }

  // Go back through the lookup table.
  double outdac = (iadc-fFlipLowAdc[range]) / fAdcDacSlope[range] + fFlipPoint[range];
  
  // Add the digital offset
  outdac += fQiePedestalDac;

  // Round to a final value.
  // FIXME: Dave simply truncates this number in the lookup table, so I will too.
  // This may change.
  int final_dac = (int) ( floor(outdac) );
  
  if(final_dac < 0) final_dac = 0; // The floor of the lookup table.

  if(final_dac > 0xffff) final_dac = 0xffff; // Final check to ensure it packs into a 16 bit word.
  MSG("DetSim",Msg::kVerbose) 
    << "QIE: " << inCharge/Munits::fC << "fC -> "
    << " True DACs " << true_dacs
    << " Range: " << range
    << " ADC: " << adc
    << " -> DAC: " << final_dac << endl;
  return final_dac;
}


void
SimQieElectronics::Print( Option_t* ) const
{
  printf("SimQieElectronics: qieDoLookup           %s\n",(fQieDoLookup)?("true"):("false"));
  printf("                   qieDacCharge          %.2f fC\n", fQieDacCharge/Munits::fC);
  printf("                   qieDacPerRangeZeroAdc %.2f DACs\n",fQieDacPerRangeZeroAdc);
  printf("                   qiePedestalWidthAdc   %.2f Range-zero ADCs\n",fQiePedestalWidthAdc);
  printf("                   qiePedestalDac        %d DAC\n",fQiePedestalDac);
  printf("                   qieSparsifyThresh     %d DAC\n",fQieSparsifyThresh);
  printf("                   doLookupNonlinearity  %s \n",fDoLookupNonlinearity?("on"):("off") );
}

void   
SimQieElectronics::Config( Registry& config )
{
  // Modify the configuration.
  //
  // Use this method to set static members with the class configuration.
  config.Get("qieDoLookup",fQieDoLookup);
  config.Get("qieDacCharge",fQieDacCharge);
  config.Get("qieDacPerRangeZeroAdc",fQieDacPerRangeZeroAdc);
  config.Get("qiePedestalWidthAdc",fQiePedestalWidthAdc);
  config.Get("qieAdcRms",fQieAdcRms);
  config.Get("qiePedestalDac",fQiePedestalDac);
  config.Get("qieSparsifyThresh",fQieSparsifyThresh);
  config.Get("doLookupNonlinearity",fDoLookupNonlinearity);

  if(fQieSparsifyThresh<fQiePedestalDac) 
    MSG("DetSim",Msg::kWarning) << "WARNING: The QIE sparsification threshold " << std::endl
                                << "         is less than the digital offset!" << std::endl
                                << "         This turns off all sparsification!" << std::endl;


  SetupTables();
}

---