SimPmt.cxx

Go to the documentation of this file.

#include "SimPmt.h"
#include "Calibrator/Calibrator.h"
#include "MessageService/MsgService.h"
#include "SimAfterpulseModel.h"
#include <fstream>

CVSID("$Id: SimPmt.cxx,v 1.31 2008/11/18 17:31:02 rhatcher Exp $");
ClassImp(SimPmt)


Bool_t   SimPmt::fsPmtDoOpticalCrosstalk = true;
Bool_t   SimPmt::fsPmtDoChargeCrosstalk  = true;
Bool_t   SimPmt::fsPmtDoNonlinearity     = true;
Bool_t   SimPmt::fsPmtDoDarkNoise        = true;
Double_t SimPmt::fsPmtScaleOpticalCrosstalk = 1.0;
Double_t SimPmt::fsPmtScaleChargeCrosstalk  = 1.0;
//Added by G. Pawloski to scale Adjacent and Diagonal separately
Double_t SimPmt::fsPmtScaleAdjacentChargeCrosstalk = 1.0;
Double_t SimPmt::fsPmtScaleDiagonalChargeCrosstalk = 1.0;
Double_t SimPmt::fsPmtScaleAdjacentOpticalCrosstalk = 1.0;
Double_t SimPmt::fsPmtScaleDiagonalOpticalCrosstalk = 1.0;
Bool_t   SimPmt::fsPmtApplyGainDrift     = true;
Bool_t   SimPmt::fsPmtHamamatsuPixelNumbering = false; // true = pixel number 1-16 instead of 0-15
Double_t SimPmt::fsVaGain;               // VA ADC->charge decalibration
Double_t SimPmt::fsQieDacCharge;         // QIE ADC->charge decalibration
Bool_t   SimPmt::fsRebuildGainMap = true;       // Rebuild the dynode chain each event, or keep it constant. Should change if gains change.

SimPmt::SimPmt ( PlexPixelSpotId tube, 
                 VldContext context, 
                 TRandom* random,
                 Int_t nPixels, Int_t nSpots 
                 )
  : fNPixels(nPixels),
    fNSpots(nSpots),
    fTube(tube),
    fContext(context),
    fRandom(random),
    fSimAfterpulseModel(0)
{
  if(random == NULL) fRandom = gRandom;
 
  SimPmt::Reset(context);
}


SimPmt::~SimPmt()
{
  BucketMap_t::iterator itr(fTimeBuckets.begin()),itrEnd(fTimeBuckets.end());
  for (; itr != itrEnd; ++itr) 
    delete itr->second; 
  fTimeBuckets.clear();

  UInt_t n= fCreatedDigiPEs.size();
  for(UInt_t i=0;i<n;i++) delete fCreatedDigiPEs[i];
}

void SimPmt::Reset(const VldContext& newContext)
{
  //
  // Clear all data in this object.
  //
  MSG("DetSim",Msg::kVerbose) << "SimPmt::Reset() " << fTube.AsString() << endl;
  fContext = newContext;
  fTotalCharge = 0;
  fDynodeTime = kPmtTime_Never;
  BucketMap_t::iterator itr(fTimeBuckets.begin()),itrEnd(fTimeBuckets.end());
  for (; itr != itrEnd; ++itr) 
    delete itr->second; 
  fTimeBuckets.clear();
  UInt_t num= fCreatedDigiPEs.size();
  for(UInt_t i=0;i<num;i++) delete fCreatedDigiPEs[i];
}



void SimPmt::AddDigiPE( const DigiPE* digipe )
{
  //
  // Add a DigiPE to the PMT.
  //

  if(!digipe) return;

  double time = digipe->GetTime();
  // First, find what bucket it goes into.
  int bucket = this->TimeToBucket(time);
  
  // Put the DigiPE in the bucket.
  // This bit of hocus-pocus does several things:
  // - makes the bucket if it doesn't exist.
  // - Translates the pixel and spot numbers to the range 1..N
  int pix =  GetPixelNumber(digipe->GetPixelSpotId());
  int spot = GetSpotNumber(digipe->GetPixelSpotId());
  GetBucket(bucket).AddDigiPE(digipe, pix, spot);

  // In addition, set the dynode time.
  GetBucket(bucket).SetDynodeTime(time);
  if(time<fDynodeTime) fDynodeTime = time;
}



Float_t SimPmt::GetPe( int pixel, int spot, int bucket ) const
{
  // 
  //  Get signal PE on pixel, spot at time bucket.
  //  Use spot = 0 to get total of all spots.
  // 
  if(spot == 0 ) {
    // Get total.
    return GetBucket(bucket).GetPixelBucket(pixel).GetTotalPE();
  }

  return GetBucket(bucket).GetPixelBucket(pixel).GetPE(spot);
}

Float_t SimPmt::GetPeXtalk( int pixel, int spot, int bucket ) const
{
  // 
  //  Get signal PE on pixel, spot at time bucket.
  //  Use spot = 0 to get total of all spots.
  // 
  if(spot == 0 ) {
    // Get total.
    return GetBucket(bucket).GetPixelBucket(pixel).GetTotalPEXtalk();
  }

  return GetBucket(bucket).GetPixelBucket(pixel).GetPEXtalk(spot);
}



Float_t SimPmt::GetCharge( int pixel, int bucket ) const
{
  // 
  //  Get charge (in couloms) on pixel, at time bucket.
  // 

  return GetBucket(bucket).GetPixelBucket(pixel).GetCharge();
}


Float_t SimPmt::GetTime( int pixel, int bucket ) const
{
  // 
  //  Get charge (in coulombs) on pixel, at time bucket.
  // 

  return GetBucket(bucket).GetPixelBucket(pixel).GetTime();
}

DigiSignal* SimPmt::CreateSignal( int pixel, int bucket ) const
{
  return GetBucket(bucket).GetPixelBucket(pixel).CreateSignal();
}


Int_t SimPmt::GetTotalHitPixels( Bool_t with_xtalk ) const
{
  int tot = 0;
  int tot_xtalk = 0;
  // Iterate over all time buckets.
  
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();
    // Iterate over pixels.
      for(Int_t pix = 1; pix <= fNPixels; pix++) {      
          SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(pix);

          if(pixtb.GetTotalPE() >0)     tot++;   
          if(pixtb.GetTotalPEXtalk() >0) tot_xtalk++;
      }
  }         
  if(with_xtalk) return tot_xtalk;
  return tot;
}

Float_t SimPmt::GetTotalPe( void ) const 
{
  
  Float_t tot = 0;
  
  // Iterate over all time buckets.
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();
    // Iterate over pixels.
      for(Int_t pix = 1; pix <= fNPixels; pix++) {      
          SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(pix);
          
          tot += pixtb.GetTotalPE();
      }
  }         

  return tot;
}

Float_t SimPmt::GetTotalCharge( void ) const 
{
  
  Float_t tot = 0;
  
  // Iterate over all time buckets.
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();
    // Iterate over pixels.
      for(Int_t pix = 1; pix <= fNPixels; pix++) {      
          SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(pix);
          
          tot += pixtb.GetCharge();
      }
  }         

  return tot;
}





Float_t SimPmt::GenPoisson( Float_t lambda, Float_t x )
{

  // For the usual compilers:
  //return  exp(-lambda + x*log(lambda)-lgamma(x+1.0));

  // For ROOT:
  return  exp(-lambda + x*log(lambda)-TMath::LnGamma(x+1.0));
}


Float_t SimPmt::RandomGenPoisson( Float_t lambda )
{

  // Sanity check
  if(lambda<=0) return 0;

  // If lambda is large enough, just use a gaussian.
  if (lambda > 88) {
      return fRandom->Gaus(0,1)*TMath::Sqrt(lambda) + lambda;
  }

  // rp is an integer, chosen at random from the poisson dist'n.  
  // rp is the nearest integer to the REAL random number we want.  
  Float_t rp = fRandom->Poisson(lambda);
  
  // Now choose a random number between rp-0.5 and rp+0.5.

  Float_t ri = TMath::Nint(rp); // Just in case it's NOT an integer.
  Float_t a = ri - 0.5;
  Float_t b = ri + 0.5;

  // Find the maximum value of f() (where f is the gen poisson f'n) in
  // this range.  If the range contains lambda, the max is lambda.
  Float_t fmax;
  if((lambda>a)&&(lambda<b)) {
    fmax = GenPoisson(lambda,lambda);
  } else {
    // Find f(a) and f(b);
    Float_t fa = GenPoisson(lambda,a);
    Float_t fb = GenPoisson(lambda,b);
    if(fa>fb) fmax = fa;
    else      fmax = fb;
  }

  // Now, pick (weighted) number between (a,b]
  while(true) {
    Float_t x = fRandom->Rndm();
    x = x + a;   
    // i.e.  x = a + x*(b-a), but (b-a) =1;
    // x is a random number (non weighted) in (a,b]

    Float_t f = GenPoisson(lambda,x);

    // r is the probability that this number is good.
    // Note we chose r in the range (0,fmax), not (0,1).
    // This keeps us from haveing to re-pick in areas
    // where the function is small.
    Float_t r = fRandom->Rndm();
    r = r*fmax;

    if(r < f) return x;
  }
}


Bool_t SimPmt::GetGainAndWidth(int pixel, int spot, Double_t& outGain, Double_t &outWidth)
{

  PlexPixelSpotId psid = GetPixelSpotId(pixel,spot);
  FloatErr adcgain, adcwidth;
  Calibrator::Instance().DecalGainAndWidth(adcgain,adcwidth,psid);

  // Set the fractional width.
  outWidth =adcwidth/adcgain;
  
  float driftedGain = adcgain;
  if(fsPmtApplyGainDrift) {
    PlexHandle plex(fContext);
    PlexStripEndId seid = plex.GetStripEndId(psid);
    if(seid.IsValid()) {
      driftedGain = Calibrator::Instance().GetDriftCorrected(adcgain,seid);
      if(driftedGain != adcgain)
        MSG("DetSim",Msg::kDebug) << "Drifted gain on " << psid.AsString() 
                                  << " " << seid.AsString() 
                                  << " by " << driftedGain/adcgain << endl;      
    }
  }
  
  // convert into gain units (i.e. unitless)
  float elecgain = 1.0/fsVaGain;
  if(GetTubeId().GetElecType()==ElecType::kQIE) elecgain = fsQieDacCharge;
  outGain = driftedGain * elecgain / Munits::e_SI;

   return true;
}


// Simulate.

void SimPmt::SimulatePmt(void)
{
  //
  // This routine, it's override, and it's helpers do the actual
  // job of simulating everything there is to simulate in the PMT.

  if(fsPmtDoOpticalCrosstalk) SimulateOpticalXtalk();     // Move single PEs around for crosstalk
  else CopyPEtoPEXtalk();  // A null operation.
  SimulateCharges();                                 // Simulate the dynode chain to get anode charge.
  if(fsPmtDoDarkNoise)        SimulateDarkNoise();        // Add some charge to some pixels by dark noise.
  if(fsPmtDoNonlinearity)     SimulateNonlinearity();     // Apply the nonlinearity
  if(fsPmtDoChargeCrosstalk)  SimulateChargeCrosstalk();  // Crosstalk some charge around.
}


void SimPmt::SimulateAfterpulsing()
{
  //
  // Creates extra PE, and adds them at delayed times.
  //

  if(fSimAfterpulseModel==0) return;

  // Iterate over all time buckets.
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();

    // Iterate over pixels
    for(Int_t injPix = 1; injPix <= fNPixels; injPix++) {      
      SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(injPix);
      
      // Iterate over all spots
      for(Int_t injSpot = 1; injSpot <= fNSpots; injSpot++ ) { 
        Float_t npe = pixtb.GetPE(injSpot);
        
        if(npe>0) {
          // Do this the easy, sloppy way: roll a number for every PE.
          SimPixelTimeBucket::PeList_t pelist = pixtb.GetDigiPE(injSpot);
          SimPixelTimeBucket::PeList_t::iterator it = pelist.begin();
          SimPixelTimeBucket::PeList_t::iterator end = pelist.end();
          for( ; it!=end; it++ ) {
            const DigiPE* sourcePe = it->second;
            // Don't afterpulse the afterpulsing...? Might be correct to do so.
            if(sourcePe->IsAfterpulse()) continue;
            
            Float_t prob = fSimAfterpulseModel->ComputeAfterpulseProb(sourcePe->GetPixelSpotId(),npe);
            if(fRandom->Uniform() < prob) {
              PlexPixelSpotId outpsid;
              Double_t outtime;
              fSimAfterpulseModel->ComputeAfterpulsePixelAndTime(
                                                                 sourcePe->GetPixelSpotId(),
                                                                 sourcePe->GetTime(),
                                                                 outpsid, outtime );
              DigiPE* newpe = new DigiPE(outtime, outpsid, DigiSignal::kAfterpulse );
              this->AddDigiPE(newpe);      // Add to simulation.
              fCreatedDigiPEs.push_back(newpe); // Add to ownership.
            }
          }
          
        }
      }
    }
  }
}

void SimPmt::SimulateOpticalXtalk()
{
  //
  // Puts extra PE into pixels they don't belong.
  //
 
  // Hack!
  // This can be uncommented to make a big dump file of the crosstalk.
  //static ofstream ofile("xtalk.dat");

  // Iterate over all time buckets.
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();

    // Iterate over the talked-to pixel
    for(Int_t xpix = 1; xpix<=fNPixels; xpix++) {  
      SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(xpix);
  
      // Iterate over talking pixels.
      for(Int_t injPix = 1; injPix <= fNPixels; injPix++) {      
        if(injPix!=xpix) {
          SimPixelTimeBucket& pixtb_inj = pmttb.GetPixelBucket(injPix);

          // Iterate over all talking spots
          for(Int_t injSpot = 1; injSpot <= fNSpots; injSpot++ ) { 
            if( pixtb_inj.GetPE(injSpot) > 0) { // There is charge in this spot.
              
              Float_t pe = GenOpticalCrosstalk(injPix, // Injected pixel
                                               injSpot,  // Injected spot
                                               xpix,     // xtalk pixel
                                               pixtb_inj.GetPE(injSpot) ); // number pe.
              
              // Debugging hack continued.
              //ofile << injPix << "\t"  << injSpot << "\t" 
              //            << xpix << "\t" << pixtb_inj.GetPE(injSpot) << "\t" << pe << endl;
              if(pe>0){
                // Choose a spot at random:
                Int_t toSpot = fRandom->Integer(fNSpots)+1;
                MoveOpticalPE(TMath::Nint(pe),
                              pixtb_inj,injSpot,
                              pixtb,    toSpot );

                pixtb.SetTruthBit(DigiSignal::kCrosstalkOptical);
              }
            }
          } // Itr over inj spot
          
        } // Itr over inj pixel
      } 
    }  // itr over xtalk pixel

    // Now we're done with all pixels in the bucket.
  } // it over buckets
  // whew!

  // Copy whatever didn't crosstalk.
  CopyPEtoPEXtalk();
}

void SimPmt::SimulateCharges()
{
  fTotalCharge = 0;
  // Iterate over all time buckets.
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();

    // Iterate over pixels.
      for(Int_t pix = 1; pix <= fNPixels; pix++) {      
          SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(pix);
          
          pixtb.SetCharge(0);
          
          for(Int_t spot = 1; spot <= fNSpots; spot++) {
            
            // This bit of magic provides the Int_trinsic
            // width of the 1 PE peak in a nice way.
            Float_t charge = GenChargeFromPE(pix,spot, pixtb.GetPEXtalk(spot));
            pixtb.AddCharge(charge);
          } // spots
          fTotalCharge += pixtb.GetCharge();
          pmttb.AddTotalCharge(pixtb.GetCharge());
      } // pixels
  } // buckets.
}


void SimPmt::SimulateDarkNoise()
{
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    Int_t ibucket = it.BucketId();
    SimPmtTimeBucket& pmttb = it.Bucket();
    
    for(Int_t pix = 1; pix <= fNPixels; pix++) {      
      SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(pix);
      
      Float_t darkcharge = GenDarkNoiseCharge( pix, GetBucketDuration(ibucket) );
      pixtb.AddCharge(darkcharge);
    }
  }
  
}


void SimPmt::SimulateNonlinearity()
{
  // 
  // Simulate nonlinearity.
  // No algorithm yet.
  //
  
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();
    
    for(Int_t pix = 1; pix <= fNPixels; pix++) {      
      SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(pix);
      
      Float_t nonlin = GenNonlinearCharge( pix, pixtb.GetCharge() );
      pixtb.SetCharge(nonlin);
      
     }
   }
}


void SimPmt::SimulateChargeCrosstalk()
{
  //
  // Simulates charge crosstalk.
  //
  // Currently just a simple fraction.

  // (Static for speed's sake.)
  static float sCharge[101]; // We have no PMT type with more than 64 pixels.
  assert(fNPixels<100);

  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();
    
    // Store old charges. 
    for(Int_t pix = 1; pix <= fNPixels; pix++) {      
      sCharge[pix] = pmttb.GetPixelBucket(pix).GetCharge();
    }

    // Do crosstalk.
    for( Int_t xpix = 1; xpix <= fNPixels; xpix++ ) {
      for( Int_t injPix = 1; injPix <=fNPixels; injPix++ ) {
        
        SimPixelTimeBucket& thePixel =  pmttb.GetPixelBucket(xpix);
        Float_t xcharge = GenElectricalCrosstalk( injPix, xpix, sCharge[injPix] );
        thePixel.AddCharge(xcharge);
        // Set truth bit for significant charge only.
        if(xcharge>1.0*Munits::fC) 
          thePixel.SetTruthBit(DigiSignal::kCrosstalk);

      }
    }
  }
}


// Actual Computation.



Float_t SimPmt::GenChargeFromPE( Int_t pixel, Int_t spot, Float_t pe )
{
  //
  // Does a single charge simulation.
  //
  
  // First version, by Nathaniel.
  // Uses the generalized poisson to get the job done.

  if(pe<=0) return 0;

  // The secondary emmission ratio is directly related to the width of the 1pe peak.
  Float_t secemm = GetPixelSecondaryEmmisionRatio(pixel,spot);

  // Choose a random number from the spectrum with a peak at the number of expected 2nary pes.
  Float_t mean2ndaryPe = pe * secemm;
  Float_t secondary_pe = RandomGenPoisson(mean2ndaryPe);

  // convert this into charge.
  return (secondary_pe / secemm)
    * GetPixelGain(pixel,spot)                   // Gain of pixel
    * 1.6e2;                                     // Convert 10^6 e -> fempto Coulombs  

}

Float_t SimPmt::GenDarkNoiseCharge( Int_t /* pixel */, 
                                    Float_t /* timeWindow */ )
{
  //
  // Dark noise simulation.
  //
  // Return the charge (in coulombs) that results from opening an integration gate of duration timeWindow (in Munits).
  // (This should return zero most of the time).
  //
  
  // No simulation yet.
  return 0; 
}


Float_t SimPmt::GenNonlinearCharge( Int_t /* pixel */,
                                    Float_t inCharge) 
{
  //
  // Nonlinearity Simulation.
  //
  // Given a charge on the anode, this routine applies
  // the nonlinearity for the PMT to give a new output charge.
  // 
  // For a completely linear response, return inCharge.
  
  return inCharge; 
}

Float_t SimPmt::GenOpticalCrosstalk( Int_t /* injPixel */,
                                     Int_t /* injSpot */,
                                     Int_t /* xtalkPixel */,
                                     Float_t /* injPE */ ) 
{
  // 
  // Optical crosstalk simulation.
  //
  // If injPE photoelectrons of light are injected into injPixel and injSpot,
  // this routine returns the number of the photoelectrons that will leak into
  // the xtalkPixel.  Note that this should probably be a random integer..
  //
  
  // Return 0 for no crosstalk.
  
  return 0; 

}

Float_t SimPmt::GenElectricalCrosstalk( Int_t /* injPixel */,                                   
                                        Int_t /* xtalkPixel */,
                                        Float_t /* injCharge */ ) 
{
  //
  // Electrical crosstalk simulation.
  //
  // If inCharge femtoCoulombs of charge are seen on the anode of injPixel, then
  // this returns the quantity of charge that will be leaked to xtalkPixel.
  // This may be a random quantity, or a fixed fraction, depending on the model.

 return 0; 
}


void SimPmt::MoveOpticalPE( UInt_t npe,
                            SimPixelTimeBucket& fromPixel,
                            Int_t               fromSpot,
                            SimPixelTimeBucket& toPixel,
                            Int_t               toSpot)
{
  //
  // Moves DigiPe from one pixel to another 
  // for use with optical crosstalk.

  // For efficiency.
  if(npe<=0) return;

  // Sanity check.
  // This can happen most likely for a single PE for which the 'mean' crosstalk is maybe 0.2%.
  // This gives a non-zero chance to get a poisson number of 2. 
  if(fromPixel.GetPE(fromSpot)<npe) {
    MSG("DetSim",Msg::kDebug) << "SimPmt::MoveOpticalPe() Optical crosstalk is bigger than injected light!" << endl;
    npe = (UInt_t)(fromPixel.GetPE(fromSpot)); // Set to max.
  }

  // Get the PE lists.
  SimPixelTimeBucket::PeList_t &fromPeList = fromPixel.GetDigiPE(fromSpot);
  SimPixelTimeBucket::PeList_t &toPeList =   toPixel.GetDigiPEXtalk(toSpot);


  // Even more sanity check.. COULD happen.
  // Happens most likely for a single PE on the pixel, which gets cross-talked
  // away to two neighboring pixels. The odds of this happening are about
  // one chance in 10,000, but there are an awful lot of pixels in an event.
  if(fromPeList.size() < npe) {
    MSG("DetSim",Msg::kDebug) << "SimPmt::MoveOpticalPe() Too much crosstalk! " 
                                << fromPeList.size() << "pe available, " 
                                << npe << " crosstalk requested." << endl;
    npe = fromPeList.size();
  }
  

  // For each moved pe
  for(UInt_t i=0; i<npe; i++) {
    // Pick a pe to move:
    Int_t whichPe = fRandom->Integer(fromPeList.size());

    // Go find it.
    SimPixelTimeBucket::PeList_t::iterator it;
    it=fromPeList.begin();
    for(int i=0;i<whichPe;i++) ++it;

    if(it!=fromPeList.end()) { // Even more sanity checks.
      toPeList.insert(*it); // Move it.
      toPixel.AddPEXtalk(toSpot,1); // Update number
      
      fromPeList.erase(it);    // Get rid of old copy: it's gone.
    }
  }
}

void SimPmt::CopyPEtoPEXtalk(void)
{
  // Iterate over all time buckets.
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();

    // Iterate over pixels.
    for(Int_t pix = 1; pix <= fNPixels; pix++) {      
      SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(pix);
      
      // Iterate over spots.
      for(Int_t spot = 1; spot<= fNSpots; spot++) {
        SimPixelTimeBucket::PeList_t &fromPeList = pixtb.GetDigiPE(spot);
        SimPixelTimeBucket::PeList_t &toPeList =   pixtb.GetDigiPEXtalk(spot);

        pixtb.AddPEXtalk(spot,(Float_t)(fromPeList.size()));
        toPeList.insert(fromPeList.begin(),fromPeList.end());   
      }
    }
  }
}

// Configuration.
void SimPmt::Config( Registry& config )
{
  Int_t tbool;
  if(config.Get("pmtDoOpticalCrosstalk",tbool)) fsPmtDoOpticalCrosstalk = tbool;
  if(config.Get("pmtDoChargeCrosstalk",tbool))  fsPmtDoChargeCrosstalk = tbool;
  if(config.Get("pmtDoNonlinearity",tbool))     fsPmtDoNonlinearity = tbool;
  if(config.Get("pmtDoDarkNoise",tbool))        fsPmtDoDarkNoise = tbool;
  if(config.Get("pmtApplyGainDrift",tbool))     fsPmtApplyGainDrift = tbool;

  config.Get("pmtScaleOpticalCrosstalk",fsPmtScaleOpticalCrosstalk);
  config.Get("pmtScaleChargeCrosstalk", fsPmtScaleChargeCrosstalk);
  //Added by G. Pawloski to scale Adjacent and Diagonal separately
  config.Get("pmtScaleAdjacentChargeCrosstalk", fsPmtScaleAdjacentChargeCrosstalk);
  config.Get("pmtScaleDiagonalChargeCrosstalk", fsPmtScaleDiagonalChargeCrosstalk);
  config.Get("pmtScaleAdjacentOpticalCrosstalk", fsPmtScaleAdjacentOpticalCrosstalk);
  config.Get("pmtScaleDiagonalOpticalCrosstalk", fsPmtScaleDiagonalOpticalCrosstalk);

  config.Get("vaGain",fsVaGain);
  config.Get("qieDacCharge",fsQieDacCharge);

  if(config.Get("pmtRebuildGainMap",tbool)) fsRebuildGainMap = tbool;

  if(config.Get("pmtHamamatsuPixelNumbering",tbool))  fsPmtHamamatsuPixelNumbering = tbool;

}

void SimPmt::PrintConfig(Option_t*)
{
  printf("SimPmt Config: pmtDoOpticalCrosstalk    %s\n",fsPmtDoOpticalCrosstalk?"true":"false");
  printf("               pmtDoChargeCrosstalk     %s\n",fsPmtDoChargeCrosstalk?"true":"false");
  printf("               pmtDoNonlinearity        %s\n",fsPmtDoNonlinearity?"true":"false");
  printf("               pmtDoDarkNoise           %s\n",fsPmtDoDarkNoise?"true":"false");
  printf("               pmtScaleOpticalCrosstalk %f\n",fsPmtScaleOpticalCrosstalk);
  printf("               pmtScaleChargeCrosstalk  %f\n",fsPmtScaleChargeCrosstalk );
  //Added by G. Pawloski to scale Adjacent and Diagonal separately
  printf("               pmtScaleAdjacentChargeCrosstalk  %f\n",fsPmtScaleAdjacentChargeCrosstalk );
  printf("               pmtScaleDiagonalChargeCrosstalk  %f\n",fsPmtScaleDiagonalChargeCrosstalk );
  printf("               pmtScaleAdjacentOpticalCrosstalk  %f\n",fsPmtScaleAdjacentOpticalCrosstalk );
  printf("               pmtScaleDiagonalOpticalCrosstalk  %f\n",fsPmtScaleDiagonalOpticalCrosstalk );
  printf("               pmtApplyGainDrift        %s\n",fsPmtApplyGainDrift?"true":"false" ); 
  printf("               pmtHamamatsuPixelNumbering %s\n",fsPmtHamamatsuPixelNumbering?"true":"false");
  printf("               vaGain                   %f ADC/fC\n",fsVaGain*Munits::fC);
  printf("               qieDacCharge             %f fC/DAC\n",fsQieDacCharge/Munits::fC);
  printf("               pmtRebuildGainMap        %s\n",fsRebuildGainMap?"true":"false");
}

---