SimPmtUTM16.cxx

Go to the documentation of this file.

#include "SimPmtUTM16.h"
#include "SimPmtMaker.h"
#include "SimPmtM16CrosstalkTable.h"
#include "SimPmtM16Crosstalk.h"
#include "MessageService/MsgService.h"
#include "MessageService/MsgFormat.h"
#include "Conventions/Munits.h"
//#include "TF1.h"

#include <cmath>
using std::pow;
using std::sqrt;

#include <cassert>
// Author: Mike Kordosky : kordosky@hep.utexas.edu
//
// Default sec values taken from original sim by J.Thomas & R.Saakyan
// This is what is in gminos.  The parameters are:
//
//  double r[12]={2.4,2.4,2.4
//            1.0,1.0,1.0,
//            1.0,1.0,1.0,
//            1.0,1.2, 2.4}; // in 100kohm units, for a total of 17.8
//  double a = 0.1472;
//  double b = 0.76;
//  double hv = 755;
//
// dynode emmission coefficients are given by
// se[i]=a*pow(hv*r[i]/rtot,b)
//
// These have been pre-calculated below:
//

const double SimPmtUTM16::gkDefaultSECValues[12] = {4.9407, 4.9407, 4.9407,
                                                    2.53997, 2.53997, 2.53997,
                                                    2.53997, 2.53997, 2.53997,
                                                    2.53997, 2.91747, 4.9407};

// the gain of the tube is the product of the emission coefficients
const double SimPmtUTM16::gkDefaultGain = SimPmtUTM16::gkDefaultSECValues[0]
*SimPmtUTM16::gkDefaultSECValues[1]
*SimPmtUTM16::gkDefaultSECValues[2]
*SimPmtUTM16::gkDefaultSECValues[3]
*SimPmtUTM16::gkDefaultSECValues[4]
*SimPmtUTM16::gkDefaultSECValues[5]
*SimPmtUTM16::gkDefaultSECValues[6]
*SimPmtUTM16::gkDefaultSECValues[7]
*SimPmtUTM16::gkDefaultSECValues[8]
*SimPmtUTM16::gkDefaultSECValues[9]
*SimPmtUTM16::gkDefaultSECValues[10]
*SimPmtUTM16::gkDefaultSECValues[11];

// global xtalk scaling parameters
double SimPmtUTM16::gOpticalXTScale=1.75;
double SimPmtUTM16::gOpticalXTOffset=0.0;
double SimPmtUTM16::gElectricXTScale=0.75;
double SimPmtUTM16::gElectricXTOffset=-0.001;

// global timing parameters
double SimPmtUTM16::gTransitTimeNS=8.8; // in ns
double SimPmtUTM16::gTransitTimeJitterNS=0.180; // in ns

// init crosstalk table to zero
SimPmtM16CrosstalkTable* SimPmtUTM16::gCrosstalkTable=0;

// Register this class with the factory.
SimPmtMakerProxy<SimPmtUTM16> gSimPmtUTM16Proxy("SimPmtUTM16");

CVSID("$Id: SimPmtUTM16.cxx,v 1.23 2005/05/13 10:57:19 tagg Exp $");

ClassImp(SimPmtUTM16)

SimPmtUTM16::SimPmtUTM16(PlexPixelSpotId tube, 
                         VldContext context,
                         TRandom* random)  :
  SimPmt(tube,context,random,16,8),
  fLastGainValidityId(-999)
{
  
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // This is a constructor

  // R1.17: Do this once at the start.
  fSECValuesBuilt = false;

  MSG("DetSim",Msg::kDebug) << "SimPmtUTM16 Constructor, Tube :" << fTube.AsString() << endl;

}

void SimPmtUTM16::Reset( const VldContext& context )
{
  // Author: N. Tagg
  //
  // Resets the PMT for a new event at a new context.
  // Ensures the correct DB is used.

  // Reset the PMT data.
  SimPmt::Reset(context);

  // Pulls new xtalk data if available, or 
  // gets it for the first time if it doesn't yet exist.
  if(gCrosstalkTable) {
    gCrosstalkTable->Reset(context);
  } else {
    gCrosstalkTable = new SimPmtM16CrosstalkTable(context);
  }

  // R1.17 Reboot the model:
  if(fsRebuildGainMap || (!fSECValuesBuilt)) 
    fSECValuesBuilt = InitSECValues();

  if(!fSECValuesBuilt){
    MSG("DetSim", Msg::kFatal)<<"Could not initialize SEC values!"<<endl;
    assert(0); // make program bomb out... is this the right thing to do?    
  }
}

void SimPmtUTM16::SimulatePmt()
{
  // Author: Mike Kordosky : kordosky@hep.utexas.edu

  if(fsPmtDoOpticalCrosstalk) SimulateOpticalXtalk();
  else CopyPEtoPEXtalk();  // A null operation.

  if(fsPmtDoDarkNoise) SimulateDarkNoise();
  SimulateCharges(); // internally uses fElecticXtalk, fNonLinearity
  
}


void SimPmtUTM16::Print(Option_t* option) const
{
  printf("  SimPmtUTM16::Print()  -- Tube: %s\n",fTube.AsString("t"));
  
  std::string opt = option;
  Bool_t spots = (opt.find("spot") != std::string::npos);;

  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    int ibucket = it.BucketId();
    
  
    printf("  Time bucket %d:\n",ibucket);

    // output:
    // pixels:  16 15 14 13   fibres:  8 7 6
    //          12 11 10  9             5 4
    //           8  7  6  5            3 2 1
    //           4  3  2  1
  
    if(spots) { 
      printf("  Photoelectrons(Npe)\n");
      for(int row=0; row < 4; row++) {
        printf("  %5.0f %5.0f %5.0f  %5.0f %5.0f %5.0f  %5.0f %5.0f %5.0f  %5.0f %5.0f %5.0f\n",
               GetPe(16-row*4,8,ibucket),  GetPe(16-row*4,7,ibucket),  GetPe(16-row*4,6,ibucket),
               GetPe(15-row*4,8,ibucket),  GetPe(15-row*4,7,ibucket),  GetPe(15-row*4,6,ibucket),
               GetPe(14-row*4,8,ibucket),  GetPe(14-row*4,7,ibucket),  GetPe(14-row*4,6,ibucket),
               GetPe(13-row*4,8,ibucket),  GetPe(13-row*4,7,ibucket),  GetPe(13-row*4,6,ibucket));

        printf("     %5.0f %5.0f        %5.0f %5.0f        %5.0f %5.0f        %5.0f %5.0f\n",
               GetPe(16-row*4,5,ibucket),  GetPe(16-row*4,4,ibucket),
               GetPe(15-row*4,5,ibucket),  GetPe(15-row*4,4,ibucket),
               GetPe(14-row*4,5,ibucket),  GetPe(14-row*4,4,ibucket),
               GetPe(13-row*4,5,ibucket),  GetPe(13-row*4,4,ibucket));
    
        printf("  %5.0f %5.0f %5.0f  %5.0f %5.0f %5.0f  %5.0f %5.0f %5.0f  %5.0f %5.0f %5.0f\n",
               GetPe(16-row*4,3,ibucket),  GetPe(16-row*4,2,ibucket),  GetPe(16-row*4,1,ibucket),
               GetPe(15-row*4,3,ibucket),  GetPe(15-row*4,2,ibucket),  GetPe(15-row*4,1,ibucket),
               GetPe(14-row*4,3,ibucket),  GetPe(14-row*4,2,ibucket),  GetPe(14-row*4,1,ibucket),
               GetPe(13-row*4,3,ibucket),  GetPe(13-row*4,2,ibucket),  GetPe(13-row*4,1,ibucket));
           
        printf("\n");
      }
    }
  
  
    printf("\n        PEs                   PE with Xtalk         Charges (fC)\n");
    for(int row=0; row < 4; row++) {
      printf("    %4d %4d %4d %4d",
             (int)GetPe(16-row*4,0,ibucket), (int)GetPe(15-row*4,0,ibucket),
             (int)GetPe(14-row*4,0,ibucket), (int)GetPe(13-row*4,0,ibucket));
    
      printf("    %4d %4d %4d %4d",
             (int)GetPeXtalk(16-row*4,0,ibucket), (int)GetPeXtalk(15-row*4,0,ibucket),
             (int)GetPeXtalk(14-row*4,0,ibucket), (int)GetPeXtalk(13-row*4,0,ibucket));

      printf("    %6.0f %6.0f %6.0f %6.0f\n",
             GetCharge(16-row*4,ibucket)/Munits::fC, GetCharge(15-row*4,ibucket)/Munits::fC,
             GetCharge(14-row*4,ibucket)/Munits::fC, GetCharge(13-row*4,ibucket)/Munits::fC);
    }
  
  }
  
  printf("\n");
}

double SimPmtUTM16::FractionalWidthToSEC(double fw)
{
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // Given a fractional width, defined as:
  // width of 1pe peak / mean of 1pe peak
  // this routine calculates the secondary emission
  // coefficient at the first 3 stages.
  // The calculation implicitly assumes that:
  // (1) We are using the minos M16 base
  // (2) Stages with the same voltage drop have the
  //     same secondary emmisson coefficient.
  //
  // For justification of the latter point see:
  // "Photomultiplier Tube", Hamamatsu Photonics, ed H. Kume, 1994
  //      Equation 3.3
  //
  //
  // A comment from PmtSim.cxx:
  // Secondary emmission ratio is derived like this:
  //  return (gain*gain)/(width*width);
  //
  // The result given above is the first order approximation.
  // It is generally good to the ~10% level.
  // For m16s with the minos base, an approximation good to
  // approximately 1%:
  // (width/gain)^2 ~= 1/sec + 1/sec^2 + 1/sec^3
  // (See NuMI-L-661, equation 12) 
  //

  // solve cubic equation for 1/sec

  const double default_sec=-1.0;

  // precondition
  // This is really a check on what emerges from the db
  // This is a rather loose criterion.  Typically
  // width/gain ~ 1/2.
  if((fw<0.01) || (fw>3.0)){
    MAXMSG("DetSim", Msg::kWarning,20)
      <<"FractionalWidthToSEC was passed the fractional width: "<<fw
      <<" which is not a reasonable value.\n"
      <<"This routine will return "<<default_sec
      <<" as the secondary emmission coefficient."<<endl;
    
    return default_sec;
  }
  
  /*
  // find solution
  // formula taken from:
  // "Mathematical Handbook of Formulas and Tables", Spiegel(1992)
  static const double Q= -2.0/9.0;
  const double R = 7.0 + 27.0*fw*fw;
  // the discriminant D = Q^3 + R^2 >0 always:
  const double D= Q*Q*Q + R*R;
  const double S=pow(R+sqrt(D),1.0/3.0);
  const double T=pow(R-sqrt(D),1.0/3.0);
  // This yields 1 real root and 2 generally complex roots.
  // The real root is:
  const double secinv= S + T - 1.0/3.0;
  // The secondary emission coefficient is:
  */

  // taken from numerical methods
  //
  const double Q = -2.0/9.0;
  const double R = (-7.0 - 27.0*fw*fw)/54.0;
  
  double secinv=1.0/default_sec;
  if( (R*R)<(Q*Q*Q)) {
    // three real roots
    MSG("DetSim", Msg::kVerbose)<<"FractionalWidthToSEC("<<fw<<") has 3 real roots!\n"
                              <<"Will use first order value: "<<1.0/(fw*fw)<<endl;
    return 1.0/(fw*fw);
  }
  else{
    const double A = pow(fabs(R) + sqrt(R*R - Q*Q*Q),1.0/3.0);
    double B=0.0; if(A!=0.0) B=Q/A;
    secinv = (A+B) - 1.0/3.0;
    MSG("DetSim", Msg::kVerbose)<<"FractionalWidthToSEC("<<fw<<") has 1 real root="
                              <<1.0/secinv<<endl;  
  }

  const double sec=1.0/secinv;
  return sec;
  
}

void SimPmtUTM16::SimulateCharges()
{

  MSG("DetSim", Msg::kVerbose)
    <<"Calling SimPmtUTM16::SimulateCharges()"<<endl;
  
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // basic structure taken from SimPmt::SimulateCharges()
  // this proceedure simulates charge xtalk in parallel
  // with the dynode amplification process
  //
  fTotalCharge = 0;
  // Iterate over all time buckets.
  MSG("DetSim", Msg::kVerbose)<<"Iterating over time buckets now."<<endl;

  int cntr=1;
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& pmttb = it.Bucket();
    MSG("DetSim", Msg::kVerbose)<<"Working on time bucket: "<<cntr<<endl;
    cntr++;
    // zero charge in all pixels of this time bucket
    // original code did the equivalent
    // not sure if it is really needed
    for(int pix=1; pix<=fNPixels; pix++) pmttb.GetPixelBucket(pix).SetCharge(0);
    // Iterate over hit pixels in this time bucket.
    MSG("DetSim", Msg::kVerbose)<<"Iterating over hit pixels."<<endl;
    for(Int_t hitpix = 1; hitpix <= fNPixels; hitpix++) {
      MSG("DetSim", Msg::kVerbose)<<"Working on pixel: "<<hitpix<<endl;
      bool oksim=false;
      for(Int_t spot = 1; spot <= fNSpots; spot++) {
        const double hitpe = pmttb.GetPixelBucket(hitpix).GetPEXtalk(spot);
        MSG("DetSim", Msg::kVerbose)<<"Working on spot: "<<spot
                                  <<" got hitpe= "<<hitpe<<endl;
        if(hitpe<=0.0) continue;
        // these collections will hold the results of the simulation
        std::vector<double> qvec(fNPixels+1, 0.0); // charge
        std::vector<double> tvec(fNPixels+1, 0.0); // transit time
        MSG("DetSim", Msg::kVerbose)<<"Calling OneHitSimDynodeChain("
                                  <<hitpix<<", "<<spot<<", "<<hitpe
                                  <<", qvec, tvec)"<<endl;      
        oksim = OneHitSimDynodeChain(hitpix, spot, hitpe, qvec, tvec);
        // add the results to the pixel time buckets
        if(oksim){
          for(int outpix=1; outpix<=fNPixels; outpix++){
            const double generated_charge = qvec[outpix-1];
            MSG("DetSim", Msg::kVerbose)<<"Generated charge="<<generated_charge
                                      <<" on pixel "<<outpix
                                      <<endl;
            pmttb.GetPixelBucket(outpix).AddCharge(generated_charge);
            // if this is not the hit pixel, then any charge is from xtalk
            if(generated_charge>0 && outpix!=hitpix){
              // flip that xtalk bit
              pmttb.GetPixelBucket(outpix).
                SetTruthBit(DigiSignal::kCrosstalk);
            }
          } 
          pmttb.AddTotalCharge(qvec[fNPixels]);
          fTotalCharge+=qvec[fNPixels];
          // what do I do with the transit time information?
        }
        else{
          MSG("DetSim", Msg::kError)<<"Error in OneHitSimDynodeChain.\nDid not place any charge into the time bucket.\nWill try to continue."<<endl;
        }
      } // spots
    } // pixels
   } // buckets.
}

bool SimPmtUTM16::OneHitSimDynodeChain(Int_t hp, Int_t hs, Float_t hpe,
                            std::vector<double>& qvec,
                            std::vector<double>& tvec)
{
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // Starting with pes on one pixel, this routine
  // simulates the secondary emission process at each dynode
  // stage.  At each stage after the first, prior to the calculation
  // of secondary electrons, this routine generates xtalk by leaking
  // some fraction of the charge in the hit pixel (at the current stage)
  // into the other pixels of the tube.  The leakage is done in a stochasitic
  // fashion with the average fraction of charge leaked taken from a database
  // table. After the simulation of charge leakage the number of secondary
  // electrons for each pixel is calculated by sampling a poisson distribution
  // with a mean given by:
  // (charge at current stage and pixel)
  // *(secondary emission coefficient at current stage and pixel)
  //
  // The charges created in this procedure are ADDED to what is already
  // in qvec.

  // NOTE:  Much of the rescaling work can probably be moved outside 
  // of this function so as to limit the number of times it must be done.
  // MAK 2003.05.12 -- now done

  // If the function was passed zero pes then quit.
  if(hpe<=0.0) return true;

  // Set some array indices
  const int hpit=hp-1;
  const int hsit=hs-1;
  MSG("DetSim", Msg::kVerbose)<<"Calculating gains."<<endl;
  // calculate the gains at each pixel
  double gains[16];
  for(int pxit=0; pxit<16; pxit++){
    gains[pxit]=1.0;
    for(int stage=0; stage<12; stage++){
      // for the hit pixel, we can use the exact gain 
      // since the spot is known
      if(pxit==hpit) gains[pxit]*=fSECValues[hpit][hsit][stage];
      // for xtalk pixels use an average gain
      else gains[pxit]*=fAverageSECValues[pxit][stage];
    }
  }

  // Now simulate dynode amplification process
  // This includes electric xtalk and pmt nonlinearity
  
  // an array to hold the charge on each pixel
  // gets updated as this routine proceeds
  unsigned int qarray[16]={0};
  
  // look up electric xtalk fractions
  // hold the fractions in xtfrac 
  MSG("DetSim", Msg::kVerbose)<<"Looking up xtalk fractions."<<endl;
  double xtfrac[16]={0.0};
  for(int pxit=0; pxit<16; pxit++){
    if(pxit==hpit) continue;
    const double got_xtfrac=RetrieveElectricXtalkFraction(hp,hs,pxit+1);
    if(got_xtfrac<0.0) xtfrac[pxit]=0.0;
    else xtfrac[pxit]=got_xtfrac;
  }

  // first stage: amplify the hit pes
  const unsigned int begin_charge=
    (unsigned int) ( fRandom->PoissonD(hpe*fSECValues[hpit][hsit][0]) );
  qarray[hpit]=begin_charge;
  MSG("DetSim", Msg::kVerbose)<<"Charge after stage 1"<<endl;
  PrintCharge(qarray);
  // loop through the last 11 stages
  // at each stage, leak some electrons for electric xtalk
  // then do dynode multiplication
  // NOTE: leakage prior to the first dynode looks like
  // what we normally call optical xtalk... that is why this
  // algorithm starts on the second stage.
  for(int stage=1; stage<12; stage++){ 
    // first leak charge for xtalk
    if(fsPmtDoChargeCrosstalk){
      for(int pxit=0; pxit<16; pxit++){
        if(pxit==hpit) continue; // pixels do not xtalk to themselves
        if(qarray[hpit]<=0.0) break; // can't work with nothing
        // create some xtalk charge
        double working_xt_fraction
          =(xtfrac[pxit]+gElectricXTOffset)*gElectricXTScale/11.0;
        // Scale it by user.
        working_xt_fraction *= fsPmtScaleChargeCrosstalk;
        qarray[pxit]+=(unsigned int)(fRandom->PoissonD(working_xt_fraction*qarray[hpit]));
      }
    }
    // amplify charges at this stage, including any charge from xtalk
    for(int pxit=0; pxit<16; pxit++){
      if(qarray[pxit]<=0) continue;// can't work with nothing
      double current_sec=0.0;
      if(pxit==hpit) current_sec=fSECValues[hpit][hsit][stage];
      else current_sec=fAverageSECValues[hpit][stage];
      // correct for non linearity on the last stage
      if(stage==11 && fsPmtDoNonlinearity){
        current_sec=
          LastStageChargeNonLinearity(qarray[pxit],gains[pxit],current_sec);
      }
      // simulate charge amplification at this stage
      qarray[pxit]=(unsigned int)(fRandom->PoissonD(qarray[pxit]*current_sec));

    }
    MSG("DetSim", Msg::kVerbose)<<"Charge after stage "<<(stage+1)<<endl;
    PrintCharge(qarray);
  }


  // dynode/anode charge fraction
  // Teststand data shows dynode/anode ~ 35%-50%
  const double dynode_fraction=0.425;
  // ADD created charges into ouput
  for(int pxit=0; pxit<16; pxit++){
    // sum charges for dynode output
    qvec[16]+=(qarray[pxit]/dynode_fraction)*Munits::e_SI;
    // convert charge from # of electrons to C
    qvec[pxit]+=qarray[pxit]*Munits::e_SI; // add charge in C
  }

  // The timing:
  // Hamamatsu says transit time 8.8 ns with a jitter of 180 ps
  // these numbers are for single pes.
  // The input charge drives the generation of charge in the tube, 
  // so simulate for it and apply the result to all pixels.

  // what is the common currency for times in DetSim?
  // I'm guessing seconds based on the above comment regarding C vs. fC.
  const double transit_time = 
    fRandom->Gaus(gTransitTimeNS, 
                  (gTransitTimeJitterNS)/sqrt(hpe))*Munits::ns;
  for(int pxit=0; pxit<17; pxit++){
    tvec[pxit]+=transit_time;
  }

  return true;
}

double SimPmtUTM16::RetrieveElectricXtalkFraction(Int_t hp, 
                                                  Int_t hs, 
                                                  Int_t xp) const
{
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // This routine should talk to the database and get the appropriate constant.
  // It can employ this object's fContext member and 
  // the arguments to this function.
  
  double result=0.0;
  if(gCrosstalkTable){
    const SimPmtM16Crosstalk* row = gCrosstalkTable->GetRow(hp,hs,xp);
    if(row){
      result = row->GetElecFrac();
      static const MsgFormat ffmt("%9.3e");
      MSG("DetSim", Msg::kVerbose)
        <<"RetrieveElectricXtalkFraction("<<hp<<", "<<hs<<", "<<xp<<") "
        <<"will return "<<ffmt(result)<<" ."<<endl;
    }
    else{
      MSG("DetSim", Msg::kError)
        <<"RetrieveElectricXtalkFraction "
        <<"could not find a db row for:\n"
        <<"hit_pixel: "<< hp
        <<", hit_spot: "<<hs
        <<", xtalk_pixel: "<<xp
        <<"."<<endl;
    }
  }
  else{
    MSG("DetSim", Msg::kFatal)
      <<"Somehow RetrieveElectricXtalkFraction was called with "
      "gCrosstalkTable==0!\nThis is a fatal error.";
      assert(0); // bomb out
  }

  return result; // a return value <= 0.0 will result in no xtalk
}

double SimPmtUTM16::RetrieveOpticalXtalkFraction(Int_t  hp, 
                                                 Int_t hs, 
                                                 Int_t xp) const
{
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // This routine should talk to the database and get the appropriate constant.
  // It can employ this object's fContext member and 
  // the arguments to this function.

  double result=0.0;
  if(gCrosstalkTable){
    const SimPmtM16Crosstalk* row = gCrosstalkTable->GetRow(hp,hs,xp);
    if(row){
      result = row->GetOptFrac();
      static const MsgFormat ffmt("%9.3e");
      MSG("DetSim", Msg::kVerbose)
        <<"RetrieveOpticalXtalkFraction("<<hp<<", "<<hs<<", "<<xp<<") "
        <<"will return "<<ffmt(result)<<" ."<<endl;
    }
    else{
      MSG("DetSim", Msg::kError)
        <<"RetrieveOpticalXtalkFraction "
        <<"could not find a db row for:\n"
        <<"hit_pixel: "<< hp
        <<", hit_spot: "<<hs
        <<", xtalk_pixel: "<<xp
        <<"."<<endl;
    }
  }
  else{
    MSG("DetSim", Msg::kFatal)
      <<"Somehow RetrieveOpticalXtalkFraction was called with "
      "gCrosstalkTable==0!\nThis is a fatal error.";
      assert(0); // bomb out
  }

  return result * fsPmtScaleOpticalCrosstalk; // a return value <= 0.0 will result in no xtalk
}


bool SimPmtUTM16::InitSECValues()
{
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // This routine fills the arrays fSECValues, fAverageSECValues
  // It should be called once when the object is made and again
  // whenever the gain and width information change for any pixel&spot
  // on this tube.
  
  MSG("DetSim",Msg::kDebug) << "InitSECValues() Tube " << fTube.AsString() << endl;

  static const double default_gain_last=gkDefaultSECValues[3]
    *gkDefaultSECValues[4]
    *gkDefaultSECValues[5]
    *gkDefaultSECValues[6]
    *gkDefaultSECValues[7]
    *gkDefaultSECValues[8]
    *gkDefaultSECValues[9]
    *gkDefaultSECValues[10]
    *gkDefaultSECValues[11];

  // rescale secs based on gains and widths taken from db
  // loop over pixels
  for(int pxit=0; pxit<16; pxit++){
    const int pixel=pxit+1;
    // loop over fiber spots on a pixel
    double ave_width=0.0;
    double ave_gain=0.0;
    for(int spit=0; spit<8; spit++){
      const int spot=spit+1;
      double pixel_gain, fractional_width;
      GetGainAndWidth(pixel,spot,pixel_gain,fractional_width);

      const double default_fractional_width=0.5;

      if(fractional_width<0.05 || fractional_width>5.0){
        // a crazy value for the spe width
        MAXMSG("DetSim", Msg::kError,50)
          <<"InitSECValues: GetGainAndWidth("<<pixel<<", "<<spot<<")"
          <<" returned fractional width " << fractional_width
          <<" which is outside the reasonable range!\n"
          <<" Will use the width "<< default_fractional_width
          <<" in further computations."<<endl;
        fractional_width=default_fractional_width;
      }
      const double sec = FractionalWidthToSEC(fractional_width);
      double width = fractional_width * pixel_gain;
      ave_width+=(width*width);
      ave_gain+=pixel_gain;
    // The following code sets the emmission coefficient of
    // the first 3 stages to sec and then rescales the 
    // default coefficients for the last 9 stages
    // so that the gain of the tube is equal to pixel_gain.
    // Since the width is largely determined by the emission
    // coefficients of the first 3 stages this procedure will result
    // in a simulation with the requested width and mean for the 1pe peak.

    // first 3 stages determine fractional width
      fSECValues[pxit][spit][0]=sec;
      fSECValues[pxit][spit][1]=sec;
      fSECValues[pxit][spit][2]=sec;

      const double pixel_gain_first=sec*sec*sec;
      // for other stages
      const double pixel_gain_last = pixel_gain/pixel_gain_first;
      const double rescale_factor = 
        pow((pixel_gain_last/default_gain_last),1.0/9.0);
      for(int stage=3; stage<12; stage++) {
        fSECValues[pxit][spit][stage] = 
          rescale_factor*gkDefaultSECValues[stage];
      }
    } // done with loop over spots

    // calculate average sec values

    // Now, ave_gain holds the sum of the gains for the 8 spots.
    // I divide by 8 to get the average.
    ave_gain/=8.0;

    // Now, ave_width holds the sum of the (fractional_width)^2 for
    // each spot. I then divide by 8 to average and take the square root.
    ave_width/=8.0;
    ave_width=sqrt(ave_width);

    // Then, the average secondary emission coeffiecient is computed from
    // ave_width, but needs to be converted back to a fractional width:
    double ave_sec=FractionalWidthToSEC(ave_width/ave_gain);

    // The first sec values of the first 3 stages are set to ave_sec.
    fAverageSECValues[pxit][0]=ave_sec;
    fAverageSECValues[pxit][1]=ave_sec;
    fAverageSECValues[pxit][2]=ave_sec;

    // The sec values of the last 9 stages are scaled to as to produce
    // an overall gain of ave_gain.
    const double ave_pixel_gain_first=ave_sec*ave_sec*ave_sec;
    const double ave_pixel_gain_last=ave_gain/ave_pixel_gain_first;
    // Figure out a rescaling factor.  The power of 1/9
    // is needed since we will end up taking this number to the 9th
    // power when computing the gain and/or simulating the dynode 
    // cascade.
    const double ave_rescale_factor = 
      pow((ave_pixel_gain_last/default_gain_last),1.0/9.0);
    for(int stage=3; stage<12; stage++) {
      fAverageSECValues[pxit][stage] = 
        ave_rescale_factor*gkDefaultSECValues[stage];
    }
    // done with rescaling business for this pixel
  }// done with loop over pixels

  // print out SEC values
  MSG("DetSim", Msg::kVerbose)<<"Printing SEC values: "<<endl;
  for(int pxit=0; pxit<16; pxit++){
    for(int spit=0; spit<8; spit++){
      MSG("DetSim", Msg::kVerbose)
        <<"Stage SECs for pixel "<<pxit+1<<", and spot "<<spit+1<<endl;
      for(int stit=0; stit<12; stit++){
      MSG("DetSim", Msg::kVerbose)
        <<"stage "<<stit+1<<" : "<<fSECValues[pxit][spit][stit]<<"   : (avg="
        <<fAverageSECValues[pxit][stit]<<")"<<endl;
      }
    }
  }

  return true;
}

double SimPmtUTM16::LastStageNonLinearity(double q, double gain, double lsec)
{
  MSG("DetSim",Msg::kError) << "YOU ARE CALLING A DEPRECATED FUNCTION!"<<endl;

  //*** DEPRECATED FUNCTION ****   Mike Kordosky -- April 25, 2005

  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // This routine models pmt non linearity as being due to a change in
  // the secondary emission of the last dynode.  The data taken from
  // the M16 teststands shows the amount of charge measured vs. the number
  // of input PEs over the range of 0-300 input pes.  The light level
  // is determined by measuring the light level in the linear 
  // region of the tube (in practice ~30pes ) and then scaling by the 
  // known transmission values of the neutral density filters used 
  // in the stand.
  //
  // Qualitatively, for light levels > 50 pes a decrease in the amount of
  // charge per input pe is seen.  This decrease is treated as being due
  // to a change in the gain of the pmt.  The data is then processed to 
  // show the fractional change in gain vs. the number of input pes.
  // That data is then parameterized as:
  // (a) Flat for pe<=50
  // (b) 3rd order polynomial for 50<pe<=300
  // (c) For pe>300 the value at pe=300 is used

  // This routine uses the parameterization above to model the change
  // in gain by a deviation of the secondary emission coefficient of the last
  // dynode stage from its nominal value.  The parametrization is defined
  // in the member function M16NonLinearity.

  // The input parameters are:
  // q = charge input to the last stage
  // gain = nominal gain of this pmt
  // lsec = nominal sec of the last stage
  
  double result=lsec;
  const double approx_npe=q/(gain/lsec);
  if(approx_npe>50.0 && approx_npe<=300.0){
    const double frac_gain_change=M16NonLinearity(approx_npe);
    result = (1.0+frac_gain_change)*lsec;
  }
  else if(approx_npe>300.0){
    const double frac_gain_change=M16NonLinearity(300.0);
    result = (1.0+frac_gain_change)*lsec;
  }

  return result;
}

double SimPmtUTM16::M16NonLinearity(double pe)
{ 
  MSG("DetSim",Msg::kError) << "YOU ARE CALLING A DEPRECATED FUNCTION!"<<endl;
  
  //***DEPRECATED FUNCTION***  Mike Kordosky -- April 25, 2005
  //
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // Returns the fractional change in gain
  // due to pmt nonlinearity.
  // Parameterization taken from m16 teststand data.
  // See SimPmtUTM16::LastStageNonLinearity
  // for more details.
  //
  // This function is to be applied for 50<pe<300.
  
  //  static const double p[4]={5.0E-2,
  //                        -1.4E-3,
  //                        5.6E-6,
  //                        7.5E-9};

  // M. Kordosky - March 22, 2005
  // sign error in 3rd degree term... I thought I fixed it long ago!!
  // thanks to Anatael for illuminating this problem
  //
  static const double p[4]={5.0E-2,
                            -1.4E-3,
                            5.6E-6,
                            -7.5E-9};
  
  return p[0] + p[1]*pe + p[2]*pe*pe + p[3]*pe*pe*pe;
}


double SimPmtUTM16::LastStageChargeNonLinearity(double q, 
                                                double /* gain */, double lsec)
{

  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // This routine models pmt non linearity as being due to a change in
  // the secondary emission of the last dynode.  The data taken from
  // the M16 teststands and CalDet LI shows the amount of charge 
  // measured vs. the number of input PEs over the range of 0-300 input pes.
  //
  // The light level is determined by measuring the light level in the linear 
  // region of the tube (in practice ~30pes ), assuming linearity there, 
  // and then scaling by the known transmission values of the 
  // neutral density filters used in the stand or the PIN measurments.
  //
  // Qualitatively, for light levels > 50 pes a decrease in the amount of
  // charge per input pe is seen.  This decrease is treated as being due
  // to a change in the gain of the pmt.  
  
  // This routine uses the parameterization above to model the change
  // in gain by a deviation of the secondary emission coefficient of the last
  // dynode stage from its nominal value.  The parametrization is defined
  // in the member function M16ChargeNonLinearity.

  // The input parameters are:
  // q = charge input to the last stage in units of # electrons
  // gain = nominal gain of this pmt (unitless)
  // lsec = nominal sec of the last stage (unitless)
  
  double result=lsec;
  const double provisional_charge=q*lsec;
  const double frac_gain_change = M16ChargeNonLinearity(provisional_charge);
  result = (1.0+frac_gain_change)*lsec;
  
  return result;
}

double SimPmtUTM16::M16ChargeNonLinearity(double nelectrons)
{ 
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // Returns the fractional change in gain
  // due to pmt nonlinearity.
  //
  // Input: the charge, in # of electron units, at final stage
  //        typically this is a provisional charge found from the 
  //        result at the second to last stage multipled by the nominal
  //        gain of the last stage
  //
  // Parameterization taken from CalDet 2003 LI
  // thanks Anatael Cabrera
  //
  // See SimPmtUTM16::LastStageChargeNonLinearity
  // for more details.
  //
  
  // default fsVaGain = 0.375/Munits::fC
  // Munits::e_SI = 1.602176462e-19 C
  const double adc_per_electron = fsVaGain*Munits::e_SI;

  // calculate provisional adcs
  // need this because parameterization given in terms of VA ADCs
  // but we haven't digitized yet!
  const double provisional_adc = nelectrons*adc_per_electron;
  
  // calculate gain dip
  // gain_dip: change of gain at last stage owing to non linearity
  //           negative # means a loss in tube gain
  //           the gain of the last stage should be modeled as
  //           gain = nominal*(1+gain_dip)
  //
  static const double p[2]={0.0166, -4.51e-6};
  
  double gain_dip = p[0] + p[1]*provisional_adc;

  const double linear_adc = 3600;//value below which PMTs are linear
  const double high_adc = 60*300;//upper limit of the parameterization ~300 PEs
  if(provisional_adc<=linear_adc) { 
    MSG("DetSim",Msg::kDebug) << " Tube in linear region: "
                              <<" provisional adc = "<<provisional_adc
                              <<"  < linear_adc = "<<linear_adc<<"\n"
                              <<" gain_dip set to zero "<<endl;
    gain_dip=0.0;
  }
  else if(provisional_adc>high_adc){
    gain_dip = p[0] + high_adc*p[1];
    MSG("DetSim",Msg::kDebug) << " Tube outside of parameterized region: "
                              <<" provisional adc = "<<provisional_adc
                              <<"  < high_adc = "<<high_adc<<"\n"
                              <<" gain_dip set to value at"
                              <<high_adc<<" ADCs : "<<gain_dip<<endl;
    
  }

  MSG("DetSim",Msg::kDebug) << "[nelectrons, provisional adc, %gain_dip] : [ " 
                            << nelectrons<<" , "
                            << provisional_adc<<" , "
                            << gain_dip*100.0 <<" ]"<<endl;


  return gain_dip; 

}


Float_t SimPmtUTM16::GenDarkNoiseCharge( Int_t /* pixel */, Float_t /* timeWindow */ )
{
  // MAK April 25, 2005 : Simulated elsewhere
  return 0; 
}

Float_t SimPmtUTM16::GenOpticalCrosstalk( Int_t hp, Int_t hs,
                                   Int_t xp, Float_t hpe ) 
{
  // Author: Mike Kordosky : kordosky@hep.utexas.edu
  //
  // This routine generates optical xtalk.  The output of the
  // routine is the number of xtalk pes generated in the xtalk pixel.
  //
  // For the configuration described by:
  // hp = hit pixel
  // hs = hit spot
  // xp = xtalk pixel
  // the database is contacted to obtain the parameter:
  // k = <xtalk_charge/hit_charge>. 
  //
  // These k values are taken from pmt test stand data and have 
  // been tuned to CalDet beam data, via the use of the global 
  // scaling parameters gOpticalXTScale, gOpticalXTOffset.
  //
  // Optical xtalk is stochastic in nature. Most of the time no charge is
  // produced. When xtalk is seen the charge corresponds to ~1pe.  
  // The average charge in the xtalk pixel is given by:
  //       average_xtalk_charge = k*hpe
  // with  k = <xtalk_charge/hit_charge>
  //
  // The number of PEs from xtalk is then generated by sampling a Poisson
  // distribution with a mean = average_xtalk_charge.

  // Can't work with nothing.
  if(hpe<=0.0) return 0.0;
  
  // Get a k value
  double k = RetrieveOpticalXtalkFraction(hp, hs, xp);
  // check k before going on
  if(k<=0.0) return 0.0; 
  
  // rescale k
  k=(k+gOpticalXTOffset)*gOpticalXTScale;
  
  // check k before going on
  if(k<=0.0) return 0.0; 
  
  // generate xtalk pes
  Int_t xpe = fRandom->Poisson(k*hpe);
  if(xpe<0) xpe=0; // some error from fRandom?
  
  return xpe;
  
}



//void SimPmtUTM16::SimulateOpticalXtalk()
//
// Puts extra PE into pixels they don't belong.
//
// MAK 2003.05.15 -- original code put all xtalk on spot 1
// fixed it to randomly generate across the tube face.
//
//
// NJT 2003.07.15 -- Put the above mod into the SimPmt class,
// remove this override.
//

void SimPmtUTM16::PrintCharge(unsigned int* qarray){
  static const MsgFormat ifmt("%10d");
  
  MSG("DetSim", Msg::kVerbose)
    <<"================================================================================\n"
    <<ifmt(qarray[3])<<" "<<ifmt(qarray[2])<<" "<<ifmt(qarray[1])<<" "<<ifmt(qarray[0])
    <<"\n"
    <<ifmt(qarray[7])<<" "<<ifmt(qarray[6])<<" "<<ifmt(qarray[5])<<" "<<ifmt(qarray[4])
    <<"\n"
    <<ifmt(qarray[11])<<" "<<ifmt(qarray[10])<<" "<<ifmt(qarray[9])<<" "<<ifmt(qarray[8])
    <<"\n"
    <<ifmt(qarray[15])<<" "<<ifmt(qarray[14])<<" "<<ifmt(qarray[13])<<" "<<ifmt(qarray[12])
    <<"\n"
    <<"================================================================================\n"

    <<endl;
}

---