SimVaTimedElectronics.cxx

Go to the documentation of this file.

#include <TMath.h>
#include "MessageService/MsgService.h"
#include "Digitization/DigiSignal.h"
#include "SimVaTimedElectronics.h"
#include "SimElecMaker.h"
#include "Calibrator/Calibrator.h"

// Register this class with the factory.
SimElecMakerProxy<SimVaTimedElectronics> gSimVaTimedElec("SimVaTimedElectronics");

CVSID("$Id: SimVaTimedElectronics.cxx,v 1.10 2007/02/07 16:08:33 rodriges Exp $");

ClassImp(SimVaTimedElectronics)

SimVaTimedElectronics::SimVaTimedElectronics(  VldContext context, TRandom* random )
  : SimVaElectronics( context, random ),
    fVaShapingTime(500.*Munits::ns)
{
}

void SimVaTimedElectronics::Config( Registry& config )
{
  // Modify the configuration.
  //
  // Use this method to set static members with the class configuration.
  config.Get("vaShapingTime",fVaShapingTime);
  SimVaElectronics::Config(config);
}


void SimVaTimedElectronics::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;
  
  // Put our PMTs into pseudo-VARCs (vmms), then deal with each varc (or vmm).
  std::map<int,vector<SimPmt*> > varcs;
  std::map<int,vector<SimPmt*> >::iterator varcItr;

  SimPmtList::iterator PmtIt;
  for(PmtIt = pmtList.begin(); PmtIt!= pmtList.end(); PmtIt++)  {
    SimPmt* pmt = PmtIt->second;
    if(pmt) {
      RawChannelId rcid = plex.GetRawChannelId(pmt->GetPixelSpotId(1));
      if ( (rcid.GetElecType()==ElecType::kVA) ) {
        if(DynodeTrigger(pmt->GetDynodeCharge())>0.) { 
          // PMT has at least 1 trigger.
          int crate = rcid.GetCrate();
          int varc  = rcid.GetVarcId();
          int vmm   = rcid.GetVmm();
          int index = crate*100 + varc*10;
          if(fVarcTriggerMode == SimVarcTriggerMode::k2of6) {
            index += vmm;
          }
          varcs[index].push_back(pmt);
        }
      }
    }
  }
  
  // Now read out each varc.
  for(varcItr=varcs.begin();varcItr!=varcs.end();varcItr++) {
    // Start with the first 
    MSG("DetSim",Msg::kDebug) << "Reading out varc " << varcItr->first << endl;
    ReadoutVarc(varcItr->second);
  }
}

void SimVaTimedElectronics::ReadoutVarc( std::vector<SimPmt*>& pmts )
{
  
  // Efficient precheck: there must be at least two pmts fired to make it go.
  if(fVarcTriggerMode != SimVarcTriggerMode::kNone) {
    if(pmts.size()<2) return; 
  };


  // Make a recording of when any dynode triggers happened.

  std::vector<double> vDynodeTimes;
  for(UInt_t i=0; i<pmts.size(); i++) {
    if(pmts[i]==0) continue; // Sanity check.
    
    // Go through all buckets in order.
    // Record dynode hits as they happen. No dead-chip time here, 
    // this is just when the ASD-Lite fires.
        
    double lastHitTime = -1e9;

    SimPmtBucketIterator bIt( *(pmts[i]) );
    for( ; !bIt.End(); bIt.Next() ) {
      Int_t bucket = bIt.BucketId();
      double bucketcharge = pmts[i]->GetDynodeCharge(bucket);
      if(DynodeTrigger(bucketcharge)) {
        double t = pmts[i]->GetDynodeTime(bucket);
        if((t-lastHitTime)>fVaChipDeadTime) {
          lastHitTime = t;
          vDynodeTimes.push_back(t);
          MSG("DetSim",Msg::kDebug) << "Dynode trigger at " << t << endl;
        }
      }
    }
  }
  
  // Now we can construct the varc trigger times.

  // Sort them.
  std::sort(vDynodeTimes.begin(),vDynodeTimes.end());

  // Vector of varc trigger times.
  std::vector<double> vVarcTrigTimes;

  if(fVarcTriggerMode != SimVarcTriggerMode::kNone) {
    for(UInt_t i=0; i<(vDynodeTimes.size() - 1); i++) {
      if((vDynodeTimes[i+1]-vDynodeTimes[i])<fVarcTriggerWindow) {
        vVarcTrigTimes.push_back(vDynodeTimes[i]);
      }
    }
  } else {
    // No trigger being asserted, so let everything sail through.
    vVarcTrigTimes = vDynodeTimes; 
  }

  // No go through the process again, looking for valid pmt triggers:
  for(UInt_t i=0; i<pmts.size(); i++) {
    if(pmts[i]==0) continue; // Sanity check.
    
    // Go through all buckets in order.
    // Do readouts as they happen, applying chip-recovery-time
    
    // Start off randomly dead, according to the formula given:
    double lastHitTime = -1e9;    

    SimPmtBucketIterator bIt( *(pmts[i]) );
    for( ; !bIt.End(); bIt.Next() ) {
      Int_t bucket = bIt.BucketId();
      double bucketcharge = pmts[i]->GetDynodeCharge(bucket);
      
      // Did it have enough charge to fire?
      if(DynodeTrigger(bucketcharge)) {
        double t = pmts[i]->GetDynodeTime(bucket);

        // Has it fired recently?
        if((t-lastHitTime)>fVaChipDeadTime) {

          // The chip was alive. 
          // Was there a valid varc trigger recently?
          bool fired = false;
          if(fVarcTriggerMode == SimVarcTriggerMode::kNone) {
            fired=true;
          }
          else {
            for(UInt_t j=0;j<vVarcTrigTimes.size();j++) {
              Double_t dt = t - vVarcTrigTimes[j];
              // if dt==0, we have the first hit that triggered the 2/36,
              // which should be read out. Use the "dt > -epsilon"
              // comparison to avoid using "==" with floats
              if ( ( dt > -1e-12 ) && ( dt < fVarcTriggerWindow ) ) {
                fired = true;
              }
            }
          }
          if(fired) {
            MSG("DetSim",Msg::kDebug) << "Chip fired. " << t/Munits::ns << endl;
            lastHitTime = t;
            ReadoutPmt(pmts[i],t);                       
          }
        }
      }
    }
  }
  
  
}



void SimVaTimedElectronics::ReadoutPmt( SimPmt* pmt, 
                                        Double_t dynodeTrigTime)
{
 MSG("DetSim",Msg::kDebug) << "SimVaTimedElectronics::ReadoutPmt " 
                           << pmt->GetTubeId().AsString() 
                           << " type " << pmt->GetType()
                           << " at time " << dynodeTrigTime
                           << endl;

  PlexHandle plex(fContext);

  // Interesting statistics:
  AddDigitsAfterFETrigger( pmt->GetTotalHitPixels(true) );
  AddAdcsAfterFETrigger( pmt->GetTotalCharge() * fVaGain );
 
  for(int ipixel = 1; ipixel<= pmt->GetNumberOfPixels(); ipixel++) {
    RawChannelId rcid = plex.GetRawChannelId(pmt->GetPixelSpotId(ipixel));
    
   int errorbits = 0;

    if(!(rcid.IsNull())) {
      //
      //  This stuff adds up all the buckets for a pixel into a single
      //  response. The charge is shaped by the VA shaping circuit.
      //

      DigiSignal* totalSignal = new DigiSignal();
      AddSignal(totalSignal);

      std::vector<Float_t> spotpe(pmt->GetNumberOfSpots()+1,0);
      
      // Find effective charge by convoluting signal with shaping curve.
      double qeff = 0;
      double qtot = 0;
      SimPmtBucketIterator it(*pmt);  
      for( ;  !it.End(); it.Next() ) { 
        int bucketId = it.BucketId(); 
        SimPmtTimeBucket& pmttb = pmt->GetBucket(bucketId);
        SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(ipixel);

        // Which spot was hit?
        for(int ispot=0;ispot<=pmt->GetNumberOfSpots();ispot++)
          spotpe[ispot]+=pixtb.GetPEXtalk(ispot);
        
        double dt = pmt->BucketToStartTime(bucketId) - dynodeTrigTime;
        
        // Skip buckets before this hit.
        if(pmt->BucketToStopTime(bucketId) <= dynodeTrigTime) continue;
        
        // Skip buckets that happen long after this hit.
        if(dt > fVaChipDeadTime) break;
        
        // Add the signal for this bucket.
        DigiSignal* signal = pixtb.CreateSignal();
        totalSignal->Merge(*signal);
        delete signal;

        double q = pixtb.GetCharge();    
        double resp = 1.0;
        if(dt>0) {
          // PHillip Litchfield tells me this is a good function form for the 
          // shaping of the VA circuit: t^2 exp(-2t/T)
          double t = dt+fVaShapingTime;
          resp = t*t*exp(-2*t/fVaShapingTime)/(fVaShapingTime*fVaShapingTime*exp(-2.0));      
        }
        qeff += q*resp;
        qtot += q;
      }
      
      MSG("DetSim",Msg::kDebug) << "Shaped readout of " << rcid.AsString() 
                                << " TrueQ = " << qtot/Munits::fC
                                << " fC.  Shaped Q = " << qeff/Munits::fC 
                                << endl; 

      // Which spot was the best one for nonlinearity?
      int bigspot = 1;
      for(int ispot=1;ispot<=pmt->GetNumberOfSpots();ispot++)
        if(spotpe[ispot] > spotpe[bigspot]) bigspot = ispot;

      PlexPixelSpotId bestPsid = pmt->GetPixelSpotId(ipixel,bigspot);

      int adc = GenSimulatedADC( qeff, rcid, bestPsid );
      int tdc = GenSimulatedTDC( dynodeTrigTime, rcid );
      
      SimDigit d(  pmt->GetPixelSpotId(ipixel),
                   rcid,
                   adc,
                   tdc,
                   totalSignal,
                   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());
      }
    }
  }

}


void
SimVaTimedElectronics::Print(Option_t* option) const
{
  SimVaElectronics::Print(option);
  printf("Timing:        fVaShapingTime        %f ns\n",fVaShapingTime/Munits::ns);
} 

---