SimPmtM64Oxford.cxx

Go to the documentation of this file.

#include "SimPmtM64Oxford.h"
#include "SimPmtM64CrosstalkTable.h"
#include "SimPmtMaker.h"
#include "SimPixelTimeBucket.h"

#include "TMath.h"

#include "MessageService/MsgService.h"
//#include "MessageService/MsgFormat.h"
#include <iostream>
#include <cassert>
using std::cout;
using std::endl;

CVSID("$Id: SimPmtM64Oxford.cxx,v 1.12 2007/01/15 20:12:15 rhatcher Exp $");


SimPmtMakerProxy<SimPmtM64Oxford> gSimPmtM64OxfordProxy("SimPmtM64Oxford");

ClassImp(SimPmtM64Oxford)


// 
// SimPmtM64Oxford
// 
// Simulates an M64 using Poission statistics at each dynode.
// It also simulates PMT non-linearity and charge smearing between QIE
// buckets. 
// Optical and electrical cross-talk simulations are inherited
// from SimPmtM64. 
// Dark noise i understand to be neglible so the (no-op) method is also 
// inherited from SimPmtM64
//

//CONSTRUCTOR

SimPmtM64Oxford::SimPmtM64Oxford(PlexPixelSpotId tube, 
                                 VldContext context,
                                 TRandom* random ): 
  SimPmtM64Full(tube,context,random)
{
  // check if multiple spots are active on the tube (probably 
  // fNSpots == 1) i'm trusting that, at the very least, fNSpots is
  //  always  non-zero & positive
  fOneSpot = true;

  fStagesBuilt = false;
      
  if (fNSpots != 1) {
    fOneSpot = false;
    MSG("DetSim",Msg::kWarning) << "This PMT has " << fNSpots 
                                << " spots! Ready for a rough ride?" 
                                << endl;
  }

  //Do some safty checks... hopefully i'll rewrite sensitive bits 
  //with vectors, then we shouldn't need this. But for now: 
  if (fNPixels > 64) {
    //were in big trouble... there pixel arrays going out of bounds
    MSG("DetSim",Msg::kFatal) << "Too many pixels! /n"
                              << "Cant continue without memory carnage!"
                              << endl;
    assert(0);  //Die in dishonor
  }
  
  
  Reset(context);

}


//STATIC DATA MEMBERS...................................................

// The stage gains are the the gain at each dynode
// Mike Kordosky called these the secondary emmision ratios in the 
// UT M16 model. currently they are calculated using the standard
// divider ratio and: 
// V(PC) = 812V
// Beta = 0.890199
// Alpha = 0.09903
// g_n = alpha * (V_n) ^ beta
Bool_t SimPmtM64Oxford::fsPmtDoChargeSmear = true;

const Double_t SimPmtM64Oxford::gkDefaultStageGains[13] = {0, 6.5397, 4.5583, 4.5583,
                                                          2.4593,  2.4593, 2.4593,
                                                          2.4593,  2.4593, 2.4593,
                                                          2.4593,  2.4593, 4.5583};

//const Double_t SimPmtM64Oxford::gkDefaultSkippingGain = 13.903;
const Double_t SimPmtM64Oxford::gkDefaultSecEmmRatio = 4.5;//5.4;// Ave over all 1pe fits 
//Double_t SimPmtM64Oxford::gkDynodeSkipRate = 0.00;//0.01 // Another rough average

const Double_t SimPmtM64Oxford::gkNLWindowBefore = 1.8056 * Munits::ns;
const Double_t SimPmtM64Oxford::gkNLWindowAfter = 0.9841 * Munits::ns;
// these define the window around a pe over which other pes can affect
// the nonlinearity.  1.8056 & 0.9841 are chosen to match (by area, not width)
// a top hat (window) to that of the 1pe pulse shape, which is aproximatly
// t^3 exp(-1.6t).   

Double_t SimPmtM64Oxford::fsDynodeSkipRate = 0.00; 
Double_t SimPmtM64Oxford::fsNLThreshold = 5;
Double_t SimPmtM64Oxford::fsNonlinearityScale = 1.5;


//MEMBER FUNCTIONS......................................................

void SimPmtM64Oxford::Reset(const VldContext& context)
{
  //Do the inherited M64 reset stuff
  SimPmtM64::Reset(context);
  
  //Reset the gains
  if((!fStagesBuilt)||fsRebuildGainMap) {
    CalStageGains();
    fStagesBuilt = true;
  }
}

//......................................................................

Bool_t SimPmtM64Oxford::CalStageGains()
{
  Bool_t success = false;

  //Calculate the default total Gain 
  Double_t DefaultTubeGain = 1;  
  for (Int_t i = 1; i <= 12; i++) {
    DefaultTubeGain *= gkDefaultStageGains[i];
  }

  
  for (Int_t pix = 1 ; pix <= fNPixels ; pix++){
    
    //Initialise the stage gains for each pixels
    for (Int_t i = 0; i <= 12; i++) {
      fStageGains[pix][i] = gkDefaultStageGains[i];
    } 
    
    //Get the gain and width from the database 
    Double_t thisPixelGain = 0;
    Double_t thisPixelWidth = 0;
    
    if ( fOneSpot ) {  
      
      success = GetGainAndWidth(pix, 1, thisPixelGain, thisPixelWidth);
      //cout << "pixel: " << pix << endl;
      //cout << "Gain set to " << thisPixelGain << endl;
      //cout << "Width set to " << thisPixelWidth << endl;;
    }  

    else{      
      //Sum over active spots
      //  as far as i know, there should only ever be one active spot.
      Double_t tempGain = 0;
      Double_t tempWidth = 0;
      Int_t nActiveSpots = 0;
      
      for (Int_t spotit = 1; spotit <= fNSpots; spotit++){
        GetGainAndWidth(pix, spotit, tempGain, tempWidth);
        if (( tempGain > 0 ) && ( tempWidth > 0 )){
          thisPixelGain += tempGain;
          thisPixelWidth += tempWidth;
          nActiveSpots++;
        } //if (spot is okay)
      }  //for (spots)
      
      thisPixelGain /= nActiveSpots;
      thisPixelWidth /= nActiveSpots;
      
    }  // >1 spots (else)

    // Truncate gains to 3 * 10^6 (should be much less)
    if (thisPixelGain > 3e6){
      MSG("DetSim",Msg::kInfo  ) << "Huge gain in DB: " << thisPixelGain 
                                 << "\tTruncating to 3e6" << endl;
      thisPixelGain = 3e6; 
    }
    //Calculate secondary emission (amplification at D1) correction
    //here we are modeling variation in the width of the 1pe peak 
    //as being entirely due to an variation in the first stage gain
  
    //first compute true (as opposed to effective) SecEmmRatio
    //See NuMi NOTE-Scint-934 for details
    
    Double_t secEmmCrtn = 1;
    Double_t effSecEmmRatio = 1/( thisPixelWidth*thisPixelWidth );   //(rho)
    
    //quick(er) method (accurate to ~0.6%, which is plenty)
    //secEmmCrtn = ( (effSecEmmRatio + 0.25) / (gkDefaultSecEmmRatio + 0.25) );

    //cout << "recal : "<< secEmmCrtn <<endl;
    
    //full, slow(er) method
    secEmmCrtn = effSecEmmRatio * ( 1. + sqrt( 1. + 1./effSecEmmRatio ) );
    secEmmCrtn /= gkDefaultSecEmmRatio * (1. + sqrt(1.+1./gkDefaultSecEmmRatio ) );
    //cout << "full recal : "<< secEmmCrtn <<endl;

    fStageGains[pix][1] *= secEmmCrtn; 

    //do inverse correction to other stages to renormalise gain

    secEmmCrtn = pow(secEmmCrtn, 1./11.); 
    
    //cout << "inv recal = "<< 1./secEmmCrtn << endl;

    for(Int_t i = 2; i <= 12; i++){
      fStageGains[pix][i] /= secEmmCrtn;
    } //for (dynodes 2->12)
   
    //calculate skipping gain 
    // using g(skip) = ( g(1)^(1/beta) + g(2)^(1/beta) )^beta
    // with beta = 0.890199
    
    const Double_t beta = 0.890199;  //Measured at test stand

      //cout << "D1 :" <<  fStageGains[pix][1] << endl;
      //cout << "D2 :" <<  fStageGains[pix][2] << endl;

    fSkippingGain[pix] = pow( fStageGains[pix][1], 1./ beta )
      + pow( fStageGains[pix][2], 1./ beta );

    //cout << "skipping gain :" << fSkippingGain[pix] << endl;

    fSkippingGain[pix] = pow( fSkippingGain[pix], beta ); 

    //cout << "skipping gain " << fSkippingGain[pix] << endl;


    //Re-calculate pixel gain

    Double_t gainRatio = 0;
    gainRatio = thisPixelGain / DefaultTubeGain ;
    //cout << "gainratio " << gainRatio << endl;
    
    gainRatio = pow(gainRatio, 1./12.);
    
    //apply correction
    
    for (Int_t i = 0; i <= 12; i++) {
      fStageGains[pix][i] *= gainRatio;
    } 

   
    //Fix the Dynode skipping gains  
    fSkippingGain[pix] =  fSkippingGain[pix] * gainRatio * gainRatio;
    
  } //for (pix)
  
  return success;
}

//......................................................................

Float_t SimPmtM64Oxford::GenChargeFromPE( Int_t pixel, 
                                          Int_t spot, Float_t pe )
{
  if (spot != 1) {
    MSG("DetSim",Msg::kError) << "Amplifiying a pe in spot "<< spot 
                              << "! This is okay... /n"
                              << "But Non-linearity will fail" 
                              << endl;
  } // spot != 1

  if (pe <= 0) return 0;  
  
  Int_t nPe = (Int_t) pe;

  if (nPe > 1200){
    //Scream like a baby
    MSG("DetSim",Msg::kInfo  ) << "BIG pulse seen: " << nPe 
                               << " photoelectons."
                               << "\tTruncating to 1200" << endl;
    nPe = 1200;
  }
  //this process is divided into two stages: 
  //the first dynodes, where individual photons count
  //the later dynodes, where we don't notice individual photons
  
  
  UInt_t neFromD2 = SimFirstDynodes(pixel, nPe);
  UInt_t neToAnode = SimLaterDynodes(pixel, neFromD2);
  
  return neToAnode * Munits::e_SI;

}

//......................................................................


UInt_t SimPmtM64Oxford::SimFirstDynodes(Int_t pixel, Int_t nPe)
{
  //bah. maybe i have to burn random numbers in a binomial to 
  //do this transparently
  //i'm pretty sure this does the same thing though 

  //Simulate 'skipping' path
  Float_t meanSkipping = nPe * fsDynodeSkipRate;
  UInt_t neThatSkip = 
    fRandom->Poisson( meanSkipping * fSkippingGain[pixel] );
  
  //simulate 'normal' path
  Float_t meanElectrons = 6.5;
  UInt_t neFromPrevDynode = nPe ;  //values for saftey

  //at dynode 1
  meanElectrons = fStageGains[pixel][1] * nPe * (1 - fsDynodeSkipRate);
  neFromPrevDynode = (UInt_t) fRandom->PoissonD(meanElectrons);

  //& at dynode 2
  meanElectrons = fStageGains[pixel][2] * neFromPrevDynode;
  neFromPrevDynode = (UInt_t) fRandom->PoissonD(meanElectrons);
                              
  //add in skipping electrons
  neFromPrevDynode += neThatSkip;

  return neFromPrevDynode;
}


//......................................................................

UInt_t SimPmtM64Oxford::SimLaterDynodes(Int_t pixel, Int_t neToD3)
{
//   Int_t switchPoint = 6;
//   //Change this if you want. It must be >= 3  and <= 12 
//   //speed = low, accuracy = high (recommend >5)


//   //At this dynode, switch from simulating each dynode individually
//   // to simulating everything in one go.
//   //It only works for Gaussian statistics!
//   //The main speed issue is with generating Poisson statisics,
//   //so bigger speed increases can be acheived by switching 
//   //at or before dynode 5. But this nessicerily gives bigger errors.
//   //There's no such thing as a free lunch!

//   //saftey net: switchPoint MUST stay within these bounds.
//   if (switchPoint < 3) switchPoint = 3;  
//   if (switchPoint > 12) switchPoint = 12;  
 
//   Float_t meanElectrons = 29.3;  //Value is just for saftey
//   UInt_t neFromPrevDynode = neToD3;
  
//   for (Int_t i = 3; i < switchPoint ; i++) {
//     meanElectrons = fStageGains[pixel][i] * (Float_t) neFromPrevDynode;
//     neFromPrevDynode = fRandom->Poisson(meanElectrons);
//    }

//   //calculate gaussian mean & width. 
//   // see NuMI note 661 for detail relating to the `extraWidth' calculation
  
//   Float_t extraWidth = 1;
//   Float_t gausGain = fStageGains[pixel][switchPoint];
  
//   for (Int_t i = (switchPoint + 1) ; i <= 12; i++ ) {
//     extraWidth = 1. + ( extraWidth * fStageGains[pixel][i] ); 
//     gausGain *= fStageGains[pixel][i];
//   }
  
//   //do gaussian amplification

//   Float_t gausMean = neFromPrevDynode * gausGain;
//   Float_t gausWidth = sqrt(gausMean * extraWidth);
  
//   neFromPrevDynode = UInt_t(fRandom->Gaus(0,1) *gausWidth + gausMean + 0.5 );
  
//..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  ..  
  //  If you want to do everything properly this does the job:

  Float_t meanElectrons = 29.3;  //Value is just for saftey
  UInt_t neFromPrevDynode = neToD3;
  for (Int_t i = 3; i <= 12; i++) {
    meanElectrons = fStageGains[pixel][i] * (Float_t) neFromPrevDynode;
    neFromPrevDynode = (UInt_t) fRandom->PoissonD(meanElectrons);
  }
  
  return neFromPrevDynode;
}

//......................................................................
void SimPmtM64Oxford::SimulateAnodeEffects()
{
  //Saftey check: 
  //if there are several spots active then there will be problems  
  Int_t activeSpot = 1; //the only spot considered in NL corrections

  if (!fOneSpot) {
    MSG("DetSim",Msg::kError  ) << "This tube has" << fNSpots 
                                << " spots active! Only PEs in spot " 
                                << activeSpot << " will be considered."
                                << endl;
  } 

  // This uses the time structure of the pe lists to calculate the
  // amount of charge active in the pmt in a window near each pe.
  // this is used to calculate a non linearity correction, which is 
  // subtracted from the pixel bucket charge.

  // The charge is then smeared into nearby buckets using a charge pulse
  // shape t^3 exp[-1.6*t]. 

  // make something to store the pe non-linearity digests
  typedef std::vector<PeNLDigest_t> AllPePixelList_t;
  
  std::vector<AllPePixelList_t> allPeList(fNPixels + 1);  
  // zeroth element is empty, so that AllPeList[n] corresponds to 
  // pixel # n, using pixels numbered 1 -> 64 ( NOT 0 ->63 )

  // Iterate over all time buckets 
  // in each time bucket: add DigiPE summaries to the lists 
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    Int_t bucketId = it.BucketId();
    SimPmtTimeBucket& pmtBucket = it.Bucket();
    
    // iterate over pixels 
    for(Int_t pix = 1 ; pix <= fNPixels ; pix++) {
      SimPixelTimeBucket& pixelBucket = pmtBucket.GetPixelBucket(pix);
        Float_t q = 0.;
        if ( pixelBucket.GetTotalPEXtalk() ) {
           q = pixelBucket.GetCharge(); 
           q /= pixelBucket.GetTotalPEXtalk(); 
        }

      //get the pe list associated with this pixel in this time bucket 
      SimPixelTimeBucket::PeList_t& peInThisBkt 
        = pixelBucket.GetDigiPEXtalk(activeSpot);  
      // spot issues... And this time, i cant iterate over spots because
      // that wouldn't be time ordered!
      
      SimPixelTimeBucket::PeList_t::iterator peItr 
        = peInThisBkt.begin();
      
      for (  ; peItr != peInThisBkt.end() ; peItr++ ) {
        
        //make a pe summary
        PeNLDigest_t peInfo(bucketId, peItr->first, q);
        //cout << peInfo.fBucketNo << endl;
  
        //store it in the relevent vector
        allPeList[pix].push_back(peInfo); 

      } //for (pes)
    } //for (pixels)
  } //for (pmt time buckets)
  
  MSG("DetSim",Msg::kVerbose) << "Constructed pe list for non-linearity"
                              << endl;


  // Now we have nice long lists of pe's spanning all time buckets
  // lets nonlinear-ise!
  
  // iterate over pixels
  for(Int_t pix = 1 ; pix <= fNPixels ; pix++) {
    
    // iterate over pes
    AllPePixelList_t::iterator currPeItr = allPeList[pix].begin();
    for (  ; currPeItr != allPeList[pix].end() ; ++currPeItr) {

      //Get the bucketId and some references.
      //Need these for Nonlinearity & charge smearing
      Int_t bucketId = currPeItr->fBucketNo;
      SimPmtTimeBucket& pmtBkt = GetBucket(bucketId);
      SimPixelTimeBucket& pixelBucket = pmtBkt.GetPixelBucket(pix);
      
      Double_t chargeRatio = 1; //if no nonlinearity

      
      if (fsPmtDoNonlinearity){ 
        //Non-linearity
        // count up how much charge is nearby
        Float_t peTime = currPeItr->fTime; 
        Float_t nearbyCharge = 0;
        AllPePixelList_t::iterator  nearbyPeItr = currPeItr;
        Float_t dt = 0;
        
        // look back
        while ( nearbyPeItr != allPeList[pix].begin() ) {
          --nearbyPeItr;
          dt = peTime - nearbyPeItr->fTime;
          if ( dt > gkNLWindowBefore ) break;
          nearbyCharge += nearbyPeItr->fCharge;
          
        } //while (past pes)
        
        //look forward, include current pe
        nearbyPeItr = currPeItr;
        
        while (nearbyPeItr != allPeList[pix].end() ) {
          dt = nearbyPeItr->fTime - peTime;
          if ( dt > gkNLWindowAfter ) break;
          nearbyCharge += nearbyPeItr->fCharge;
          ++nearbyPeItr;
        }  //while (future pes)
        
        
        // calculate NL correction based on this amount
        // newCharge is the nonlinearised charge for all the pes in
        // the window, we want it for just one.
        if (nearbyCharge){
          // because it is possable have pes but no charge
          // (no seconday emission at some stage)
          chargeRatio = GenNonlinearCharge( pix, nearbyCharge );
          chargeRatio /= nearbyCharge;
        }       
                
        // use this correction to reduce the charge in this bucket
        Float_t deltaq = (chargeRatio - 1) * currPeItr->fCharge;
        pixelBucket.AddCharge(deltaq);
      
      } // if (fsDoNonlinearity) 



      if (fsPmtDoChargeSmear){
        //Smear the charge about
        
        // "This smushes charge into the next and previous QIE buckets.  
        // It's assumed that digiPE occours at the peak of the pulse, 
        // which has a shape of t^3 exp(-t*1.6). 
        // Some of the charge gets put into the next and previous buckets.
        // Note that the truth info is lost: you will now get digits with 
        // no signal info... but that's the way it goes, I'm afraid."
        
        Float_t startTime = BucketToStartTime(bucketId);
        Float_t stopTime = BucketToStopTime(bucketId);
        
        SimPmtTimeBucket& nextPmtBucket = GetBucket(bucketId + 1);
        SimPixelTimeBucket& nextPixBucket = nextPmtBucket.GetPixelBucket(pix);
        SimPmtTimeBucket& prevPmtBucket = GetBucket(bucketId - 1); 
        SimPixelTimeBucket& prevPixBucket = prevPmtBucket.GetPixelBucket(pix);
        
        //How much into the adjacent buckets?
        Float_t qNext = currPeItr->fCharge * chargeRatio;
        qNext *= NextBucketFraction ( stopTime - currPeItr->fTime );
        Float_t qPrev = currPeItr->fCharge * chargeRatio;
        qPrev *= PrevBucketFraction ( currPeItr->fTime - startTime );
        
        //cout << "Time from boundary = " << (currPeItr->fTime - startTime) * 1e9 <<endl;
        
        //Add charge to next bucket, and set truth bits if there is a 
        //non-trivial amount.
        nextPixBucket.AddCharge(qNext);
        if ( qNext > 1.0 * Munits::fC ) {
          nextPixBucket.SetTruthBit(DigiSignal::kLeakFromPrevBucket);
        }
        
        //Add charge to previous bucket, and set truth bits if there is a 
        //non-trivial amount.
        prevPixBucket.AddCharge(qPrev);
        if ( qPrev > 1.0 * Munits::fC ) {
          prevPixBucket.SetTruthBit(DigiSignal::kLeakFromNextBucket);
        }
        
        //remove the same amount of charge from the current bucket
        pixelBucket.AddCharge( -(qPrev + qNext) ); 
      
      } // if (fsDoChargeSmear)
    } // for (pes)
  } // for (pixels)
  
  MSG("Detsim",Msg::kVerbose) << "Nonlinearity and Charge smearing done"
                              << endl;
    
}
//......................................................................

Float_t SimPmtM64Oxford::GenNonlinearCharge( Int_t /*pixel*/,
                                             Float_t linCharge) 
{
  //Debug
  //if (fsNonlinearityScale == 0) {cout << "fsnls: exactly 0" << endl;}
  //else {cout << "fsnls: " << fsNonlinearityScale << endl;} 

  //Note on Tuning parameters: threshold(fsNLThreshold) and dGdQ
 
  //Threshold in terms of (approx) number of pes in window.
  //For pulses similar in shape to the oxford test stand, and a ~3ns window
  //the non-linearity `turns on' for pulses sizes (in pe) of ~7 times
  //this threshold.
  
  //dGdQ is the (Default) rate at which the gain drops off at threshold
  //increasing fsNonlinearityScale -> NL turns on more quickly

  Float_t scaleNL = 1.;
  const Float_t qInPe = linCharge / (1e6 * Munits::e_SI);  

  const Float_t dGdQ = 3.9e-3; //Best tune from Test Stand  

  //Exponential drop off model:
  if (qInPe > fsNLThreshold) {
    const Float_t alpha = 2 * dGdQ * fsNonlinearityScale;  
    //makes everything  easier to read. 
    
    scaleNL = 1 - TMath::Exp( -alpha * (qInPe - fsNLThreshold) );
    scaleNL /= alpha * ( qInPe - fsNLThreshold );  

  }  // if (above threshold)

return linCharge * scaleNL;
  
}
//......................................................................

void SimPmtM64Oxford::SimulateCharges()
{
  fTotalCharge =0;
  
  //cout << "PMT total pe " << GetTotalPe() << endl; 

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

    //Iterate over pixels
    for(Int_t pix = 1 ; pix <= fNPixels ; pix++) {
      SimPixelTimeBucket& pixelBucket = pmtBucket.GetPixelBucket(pix);
      Float_t nPe = pixelBucket.GetTotalPEXtalk();
      
      //Generate a charge in the pixel bucket
      //iterate over spots (...)
      for (Int_t spot = 1 ; spot <= fNSpots ; spot++){
        Float_t q = GenChargeFromPE( pix, spot, nPe );
        //cout << "Anode charge = " << q << endl;
        pixelBucket.AddCharge(q);
        fTotalCharge += q;
      } //for (spots)
    } //for (pixels)
  } // for (pmt time buckets)
}

//----------------------------------------------------------------------

void SimPmtM64Oxford::Config( Registry& config )
{
  double dtmp;
  int itmp;
  if(config.Get("pmtDoChargeSmear",itmp)) fsPmtDoChargeSmear = itmp;
  if(config.Get("pmtM64DynodeSkipRate",dtmp)) fsDynodeSkipRate = dtmp;
  if(config.Get("pmtM64NLThreshold",dtmp)) fsNLThreshold = dtmp;
  if(config.Get("pmtM64NonlinearityScale",dtmp)) fsNonlinearityScale =dtmp;

  SimPmtM64Full::Config(config);

  // A non linearity scale of zero is equivalent to no nonlinearity,
  // so turn it off now because it would also cause a FPE.
  // this is quicker than checking each time we go past.  
  if (0 == fsNonlinearityScale) fsPmtDoNonlinearity = 0;
}

//----------------------------------------------------------------------

void SimPmtM64Oxford::SimulatePmt()
{
  //
  // Reimplemented from SimPmt. Nonlinearity and Charge smearing are both done
  // in a new function SimulateAnodeEffects.
  //
  //cout << "*************************************************************" << endl;
  // cout << "        Start simulation." << endl;
  //cout << "*************************************************************" << endl;
  //Print();
  
  if(fsPmtDoOpticalCrosstalk) SimulateOpticalXtalk();     // Move single PEs around for crosstalk
  else CopyPEtoPEXtalk();  // A null operation.
  
  //cout << "*************************************************************" << endl;
  // cout << "        After Optical." << endl;
  //cout << "*************************************************************" << endl;
  //Print();
  
  SimulateCharges();                                 // Simulate the dynode chain to get anode charge.
  
  //cout << "*************************************************************" << endl;
  //cout << "        After Charges." << endl;
  //cout << "*************************************************************" << endl;
  //Print();
  
  if(fsPmtDoDarkNoise)        SimulateDarkNoise();        // Add some charge to some pixels by dark noise.
  if(fsPmtDoChargeSmear ||
     fsPmtDoNonlinearity   )     SimulateAnodeEffects();     // Apply the nonlinearity
  if(fsPmtDoChargeCrosstalk)  SimulateChargeCrosstalk();  // Crosstalk some charge around.
  //cout << "*************************************************************" << endl;
  //cout << "        Final." << endl;
  //cout << "*************************************************************" << endl;
  //Print();
 }

---