SimPmtM16.cxx

Go to the documentation of this file.

#include "SimPmtM16.h"
#include "MessageService/MsgService.h"
#include "SimPmtMaker.h"
#include <cmath>

// Register this class with the factory.
SimPmtMakerProxy<SimPmtM16> gSimPmtM16Proxy("SimPmtM16");

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

ClassImp(SimPmtM16)

SimPmtM16::SimPmtM16(PlexPixelSpotId tube, 
                     VldContext context,
                     TRandom* random)  :
  SimPmt(tube,context,random,16,8) 
{
}


void SimPmtM16::Print(Option_t* /*option*/) const
{
  printf("  SimPmtM16::Print()  -- Tube: %s\n",fTube.AsString("t"));

 // Iterate over all time buckets.
  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
  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("  Charges (fC)\n");
  for(int row=0; row < 4; row++) {
    // printf("    %7.1f %7.1f %7.1f %7.1f\n",
    printf("    %f \t %f \t %f \t %f\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");
  }
}


Float_t SimPmtM16::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.
  Double_t gain, width;
  GetGainAndWidth(pixel,spot,gain,width);
  
  // Secondary emmission ratio is (mean*mean)/(rms*rms). 'width' is fractional width, so:
  Float_t secemm = 1/(width*width);

  // 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. 
  Float_t charge = (secondary_pe / secemm)       // Charge in PE-units
    * gain
    * Munits::e_SI;               // e to charge
  

  MSG("DetSim",Msg::kVerbose) << "SimPmtM16: pe to charge: " 
                              << pe << "\t"
                              << secemm << "\t"
                              << secondary_pe << "\t"
                              << gain << "\t"
                              << charge << endl;
  return charge;
}

Float_t SimPmtM16::GenDarkNoiseCharge( Int_t /* pixel */, Float_t /* timeWindow */ )
{
  //
  // Dark noise simulation.
  //
  // Return the charge (in fC) 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 SimPmtM16::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; 
}


// The following is pretty much copied-and-pasted out of
// m16_xt.F, the xtalk simulation put in by Mike Kordosky
// in GMINOS.

const float m16_opt_bc[] = { 0, 4.1E-4, 30.8E-4, 7.1E-4,
                     15.7E-4,         43.8E-4,
                     4.4E-4, 44.4E-4, 8.6E-4,
                     0.0,     0.0,    0.0};
const float m16_ul_rprob[] = { 0, 1.0351, 1.40591, 2.51808,
                       1.0000, 1.3926,
                       0.939646, 1.05889, 1.33843};
const float m16_u_rprob[] = { 0, 1.92966, 2.08835, 1.93475,
                      1.0000, 1.01693,
                      0.581402, 0.58964, 0.548125};
const float m16_ur_rprob[] = { 0, 1.77202, 1.10947, 0.881364,
                       1.0000, 0.869394,
                       1.02619, 0.902475, 0.821414};
const float m16_l_rprob[] = { 0, 0.742614, 1.46936, 4.88177,
                      1.0000, 2.72021,
                      0.80697, 1.495, 4.89158};
const float m16_r_rprob[] = { 0, 1.69269, 0.626856, 0.372045,
                      1.0000, 0.489508,
                      1.68536, 0.642841, 0.391364};
const float m16_dl_rprob[] = { 0, 0.841386, 1.04104, 1.71787,
                       1.0000, 1.79258,
                       1.28237, 1.79684, 4.16197};
const float m16_d_rprob[] = { 0, 0.450814, 0.469413, 0.428049,
                      1.0000, 1.00814,
                      2.0646, 2.23663, 2.11148};
const float m16_dr_rprob[] = { 0, 0.964592, 0.65102, 0.561097,
                       1.0000, 0.815038,
                       2.25732, 1.13462, 0.816332};
const float m16_wd_rprob[] = { 0, 1.0, 1.0, 1.0,
                       1.0, 1.0,
                       1.0, 1.0, 1.0};
const float m16_pd_rprob[] = { 0, 1.0, 1.0, 1.0,
                       1.0, 1.0,
                       1.0, 1.0, 1.0};
const float m16_o_rprob[] = { 0, 1.0, 1.0, 1.0,
                      1.0, 1.0,
                      1.0, 1.0, 1.0};


Float_t SimPmtM16::GenOpticalCrosstalk( Int_t injPixel, Int_t injSpot,
                                   Int_t xPixel, 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.
    //   Pixels:
  //   4  3  2  1
  //   8  7  6  5
  //   12 11 10 9
  //   16 15 14 13

  //     8 border pixels:
  //     1 2 3
  //     4 X 5   (X = hit pixel)
  //     6 7 8
  //
  //     All others:
  //     9 = Far neighbor going with dynode structure
  //     10 = Far neighbor going perp to dynode structure
  //     11 = Far diagonal neighbor
  //     0 = hit_pixel.eq.xt_pixel ! an error
  if(injPixel==xPixel) return 0;
  
  Float_t coeff;
  Float_t rprob;

  if(xPixel==(injPixel-4)) { // 2
    coeff = m16_opt_bc[2];
    rprob = m16_u_rprob[injSpot];
  }
  else if(xPixel==(injPixel+4)) { // 7
    coeff = m16_opt_bc[7];
    rprob = m16_d_rprob[injSpot];
  }
  else if((xPixel==(injPixel+1))&&(injPixel%4!=0)) { // 4
    coeff = m16_opt_bc[4];
    rprob = m16_l_rprob[injSpot];
  }
  else if((xPixel==(injPixel-1))&&((injPixel-1)%4!=0)) { // 5
    coeff = m16_opt_bc[5];
    rprob = m16_r_rprob[injSpot];
  }
  else if((injPixel%4!=0)&&(xPixel==injPixel-3)) { //1
    coeff = m16_opt_bc[1];
    rprob = m16_ul_rprob[injSpot];
  }
  else if(((injPixel-1)%4!=0)&&(xPixel==injPixel-5)) {  // 3
    coeff = m16_opt_bc[3];
    rprob = m16_ur_rprob[injSpot];
  }
  else if((injPixel%4!=0)&&(xPixel==injPixel+5)) { // 6
    coeff = m16_opt_bc[6];
    rprob = m16_dl_rprob[injSpot];
  }
  else if(((injPixel-1)%4!=0)&&(xPixel==injPixel+3)) {  // 8
    coeff = m16_opt_bc[8];
    rprob = m16_dr_rprob[injSpot];
  }
  else if( abs(injPixel-xPixel)%4 == 0 ) { // 9
    coeff = m16_opt_bc[9];
    rprob = m16_wd_rprob[injSpot];
  } 
  else if(((injPixel-1)%4) == ((xPixel-1)%4)) { // 10
    coeff = m16_opt_bc[10];
    rprob = m16_pd_rprob[injSpot];
  }
  else{ // 11
    coeff = m16_opt_bc[11];
    rprob = m16_o_rprob[injSpot];
  }

  // correct poisson coefficient for pixel spot
  // dependence
  const Float_t ref_pe = 35.0;
  const Float_t m16_opt_xt_scale = 2.0;

  Float_t cor1_coeff = -(1.0/ref_pe)*
    log(1.0-rprob*(1.0-exp(-coeff*ref_pe)));
  // correct poisson coefficient for overall scale
  Float_t cor2_coeff = -(1.0/ref_pe)*
    log(1.0-m16_opt_xt_scale*(1.0-exp(-cor1_coeff*ref_pe)));
  // calculate xt probability

  // HACK HACK HACK
  // This adds just a wee little bit of crosstalk to any
  // pair of pixel/spots.  Total addition is ~1 percent total
  // crosstalk, but across the whole face of the pmt.
  // cor2_coeff += 0.0005;
  
  // Scale it.
  cor2_coeff *= fsPmtScaleOpticalCrosstalk;

  return (Float_t)fRandom->Poisson(cor2_coeff * injPE);
}

Float_t SimPmtM16::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; 
}

---