SimPmtUTM16 Class Reference

#include <SimPmtUTM16.h>

Inheritance diagram for SimPmtUTM16:
SimPmt SimPmtM16UTTimed

List of all members.

Public Member Functions

 SimPmtUTM16 (PlexPixelSpotId tube, VldContext context, TRandom *random=NULL)
virtual ~SimPmtUTM16 ()
virtual void Reset (const VldContext &context)
virtual Pmt_t GetType () const
virtual void Print (Option_t *option="") const

Protected Member Functions

virtual Float_t GetDynodeGain () const
virtual Float_t GenDarkNoiseCharge (Int_t pixel, Float_t timeWindow)
virtual Float_t GenOpticalCrosstalk (Int_t injPixel, Int_t injSpot, Int_t xtalkPixel, Float_t injPE)
virtual void SimulateCharges ()
virtual void SimulatePmt ()

Private Member Functions

bool OneHitSimDynodeChain (Int_t hp, Int_t hs, Float_t hpe, std::vector< double > &qvec, std::vector< double > &tvec)
double ScaleDiagonalAdjacentChargeCrosstalk (Int_t pxit, Int_t hpit) const
double ScaleDiagonalAdjacentOpticalCrosstalk (Int_t pxit, Int_t hpit) const
double RetrieveElectricXtalkFraction (Int_t hp, Int_t hs, Int_t xp) const
double RetrieveOpticalXtalkFraction (Int_t hp, Int_t hs, Int_t xp) const
void PrintCharge (unsigned int *q)
bool InitSECValues ()

Static Private Member Functions

static double FractionalWidthToSEC (double fw)
static double LastStageNonLinearity (double q, double gain, double lastsec)
static double M16NonLinearity (double pe)
static double LastStageChargeNonLinearity (double q, double gain, double lastsec)
static double M16ChargeNonLinearity (double pe)

Private Attributes

Int_t fLastGainValidityId
bool fSECValuesBuilt
double fSECValues [16][8][12]
double fAverageSECValues [16][12]
std::set< int > adjacentPixels [16]
std::set< int > diagonalPixels [16]

Static Private Attributes

static const double gkDefaultSECValues [12]
static const double gkDefaultGain
static double gOpticalXTScale = 1.75
static double gOpticalXTOffset = 0.0
static double gElectricXTScale = 0.75
static double gElectricXTOffset = -0.001
static double gTransitTimeNS = 8.8
static double gTransitTimeJitterNS = 0.180
static SimPmtM16CrosstalkTablegCrosstalkTable = 0

Detailed Description

Definition at line 46 of file SimPmtUTM16.h.


Constructor & Destructor Documentation

SimPmtUTM16::SimPmtUTM16 ( PlexPixelSpotId  tube,
VldContext  context,
TRandom *  random = NULL 
)

Definition at line 73 of file SimPmtUTM16.cxx.

References Msg::kDebug, and MSG.

00075                                            :
00076   SimPmt(tube,context,random,16,8),
00077   fLastGainValidityId(-999)
00078 {
00079   
00080   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00081   //
00082   // This is a constructor
00083 
00084   // R1.17: Do this once at the start.
00085   fSECValuesBuilt = false;
00086 
00087   //Below code added by G. Pawloski to scale diagonal and adjacent Xtalk separately 
00088   // Lookup tables of adjacent and diagonal pixels
00089   // There might be a smarter way of doing this but using a lookup table is less prone to logic errors
00090   //   Pixels:
00091   //   3  2  1  0
00092   //   7  6  5  4
00093   //   11 10 9  8
00094   //   15 14 13 12
00095   adjacentPixels[0].insert(1);
00096   adjacentPixels[0].insert(4);
00097   adjacentPixels[1].insert(0);
00098   adjacentPixels[1].insert(2);
00099   adjacentPixels[1].insert(5);
00100   adjacentPixels[2].insert(1);
00101   adjacentPixels[2].insert(3);
00102   adjacentPixels[2].insert(6);
00103   adjacentPixels[3].insert(2);
00104   adjacentPixels[3].insert(7);
00105   adjacentPixels[4].insert(0);
00106   adjacentPixels[4].insert(5);
00107   adjacentPixels[4].insert(8);
00108   adjacentPixels[5].insert(1);
00109   adjacentPixels[5].insert(4);
00110   adjacentPixels[5].insert(6);
00111   adjacentPixels[5].insert(9);
00112   adjacentPixels[6].insert(2);
00113   adjacentPixels[6].insert(5);
00114   adjacentPixels[6].insert(7);
00115   adjacentPixels[6].insert(10);
00116   adjacentPixels[7].insert(3);
00117   adjacentPixels[7].insert(6);
00118   adjacentPixels[7].insert(11);
00119   adjacentPixels[8].insert(4);
00120   adjacentPixels[8].insert(9);
00121   adjacentPixels[8].insert(12);
00122   adjacentPixels[9].insert(5);
00123   adjacentPixels[9].insert(8);
00124   adjacentPixels[9].insert(10);
00125   adjacentPixels[9].insert(13);
00126   adjacentPixels[10].insert(6);
00127   adjacentPixels[10].insert(9);
00128   adjacentPixels[10].insert(11);
00129   adjacentPixels[10].insert(14);
00130   adjacentPixels[11].insert(7);
00131   adjacentPixels[11].insert(10);
00132   adjacentPixels[11].insert(15);
00133   adjacentPixels[12].insert(8);
00134   adjacentPixels[12].insert(13);
00135   adjacentPixels[13].insert(9);
00136   adjacentPixels[13].insert(12);
00137   adjacentPixels[13].insert(14);
00138   adjacentPixels[14].insert(10);
00139   adjacentPixels[14].insert(13);
00140   adjacentPixels[14].insert(15);
00141   adjacentPixels[15].insert(11);
00142   adjacentPixels[15].insert(14);
00143 
00144   diagonalPixels[0].insert(5);
00145   diagonalPixels[1].insert(4);
00146   diagonalPixels[1].insert(6);
00147   diagonalPixels[2].insert(5);
00148   diagonalPixels[2].insert(7);
00149   diagonalPixels[3].insert(6);
00150   diagonalPixels[4].insert(1);
00151   diagonalPixels[4].insert(9);
00152   diagonalPixels[5].insert(0);
00153   diagonalPixels[5].insert(2);
00154   diagonalPixels[5].insert(8);
00155   diagonalPixels[5].insert(10);
00156   diagonalPixels[6].insert(1);
00157   diagonalPixels[6].insert(3);
00158   diagonalPixels[6].insert(9);
00159   diagonalPixels[6].insert(11);
00160   diagonalPixels[7].insert(2);
00161   diagonalPixels[7].insert(10);
00162   diagonalPixels[8].insert(5);
00163   diagonalPixels[8].insert(13);
00164   diagonalPixels[9].insert(4);
00165   diagonalPixels[9].insert(6);
00166   diagonalPixels[9].insert(12);
00167   diagonalPixels[9].insert(14);
00168   diagonalPixels[10].insert(5);
00169   diagonalPixels[10].insert(7);
00170   diagonalPixels[10].insert(13);
00171   diagonalPixels[10].insert(15);
00172   diagonalPixels[11].insert(6);
00173   diagonalPixels[11].insert(14);
00174   diagonalPixels[12].insert(9);
00175   diagonalPixels[13].insert(8);
00176   diagonalPixels[13].insert(10);
00177   diagonalPixels[14].insert(9);
00178   diagonalPixels[14].insert(11);
00179   diagonalPixels[15].insert(10);
00180 
00181   MSG("DetSim",Msg::kDebug) << "SimPmtUTM16 Constructor, Tube :" << fTube.AsString() << endl;
00182 
00183 }

virtual SimPmtUTM16::~SimPmtUTM16 (  )  [inline, virtual]

Definition at line 53 of file SimPmtUTM16.h.

00053 {};


Member Function Documentation

double SimPmtUTM16::FractionalWidthToSEC ( double  fw  )  [static, private]

Definition at line 291 of file SimPmtUTM16.cxx.

References Msg::kVerbose, Msg::kWarning, MAXMSG, MSG, and MuELoss::R.

Referenced by InitSECValues().

00292 {
00293   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00294   //
00295   // Given a fractional width, defined as:
00296   // width of 1pe peak / mean of 1pe peak
00297   // this routine calculates the secondary emission
00298   // coefficient at the first 3 stages.
00299   // The calculation implicitly assumes that:
00300   // (1) We are using the minos M16 base
00301   // (2) Stages with the same voltage drop have the
00302   //     same secondary emmisson coefficient.
00303   //
00304   // For justification of the latter point see:
00305   // "Photomultiplier Tube", Hamamatsu Photonics, ed H. Kume, 1994
00306   //      Equation 3.3
00307   //
00308   //
00309   // A comment from PmtSim.cxx:
00310   // Secondary emmission ratio is derived like this:
00311   //  return (gain*gain)/(width*width);
00312   //
00313   // The result given above is the first order approximation.
00314   // It is generally good to the ~10% level.
00315   // For m16s with the minos base, an approximation good to
00316   // approximately 1%:
00317   // (width/gain)^2 ~= 1/sec + 1/sec^2 + 1/sec^3
00318   // (See NuMI-L-661, equation 12) 
00319   //
00320 
00321   // solve cubic equation for 1/sec
00322 
00323   const double default_sec=-1.0;
00324 
00325   // precondition
00326   // This is really a check on what emerges from the db
00327   // This is a rather loose criterion.  Typically
00328   // width/gain ~ 1/2.
00329   if((fw<0.01) || (fw>3.0)){
00330     MAXMSG("DetSim", Msg::kWarning,20)
00331       <<"FractionalWidthToSEC was passed the fractional width: "<<fw
00332       <<" which is not a reasonable value.\n"
00333       <<"This routine will return "<<default_sec
00334       <<" as the secondary emmission coefficient."<<endl;
00335     
00336     return default_sec;
00337   }
00338   
00339   /*
00340   // find solution
00341   // formula taken from:
00342   // "Mathematical Handbook of Formulas and Tables", Spiegel(1992)
00343   static const double Q= -2.0/9.0;
00344   const double R = 7.0 + 27.0*fw*fw;
00345   // the discriminant D = Q^3 + R^2 >0 always:
00346   const double D= Q*Q*Q + R*R;
00347   const double S=pow(R+sqrt(D),1.0/3.0);
00348   const double T=pow(R-sqrt(D),1.0/3.0);
00349   // This yields 1 real root and 2 generally complex roots.
00350   // The real root is:
00351   const double secinv= S + T - 1.0/3.0;
00352   // The secondary emission coefficient is:
00353   */
00354 
00355   // taken from numerical methods
00356   //
00357   const double Q = -2.0/9.0;
00358   const double R = (-7.0 - 27.0*fw*fw)/54.0;
00359   
00360   double secinv=1.0/default_sec;
00361   if( (R*R)<(Q*Q*Q)) {
00362     // three real roots
00363     MSG("DetSim", Msg::kVerbose)<<"FractionalWidthToSEC("<<fw<<") has 3 real roots!\n"
00364                               <<"Will use first order value: "<<1.0/(fw*fw)<<endl;
00365     return 1.0/(fw*fw);
00366   }
00367   else{
00368     const double A = pow(fabs(R) + sqrt(R*R - Q*Q*Q),1.0/3.0);
00369     double B=0.0; if(A!=0.0) B=Q/A;
00370     secinv = (A+B) - 1.0/3.0;
00371     MSG("DetSim", Msg::kVerbose)<<"FractionalWidthToSEC("<<fw<<") has 1 real root="
00372                               <<1.0/secinv<<endl;  
00373   }
00374 
00375   const double sec=1.0/secinv;
00376   return sec;
00377   
00378 }

Float_t SimPmtUTM16::GenDarkNoiseCharge ( Int_t  pixel,
Float_t  timeWindow 
) [protected, virtual]

Reimplemented from SimPmt.

Definition at line 1087 of file SimPmtUTM16.cxx.

01088 {
01089   // MAK April 25, 2005 : Simulated elsewhere
01090   return 0; 
01091 }

Float_t SimPmtUTM16::GenOpticalCrosstalk ( Int_t  injPixel,
Int_t  injSpot,
Int_t  xtalkPixel,
Float_t  injPE 
) [protected, virtual]

Reimplemented from SimPmt.

Definition at line 1093 of file SimPmtUTM16.cxx.

References SimPmt::fRandom, gOpticalXTOffset, gOpticalXTScale, and RetrieveOpticalXtalkFraction().

01095 {
01096   // Author: Mike Kordosky : kordosky@hep.utexas.edu
01097   //
01098   // This routine generates optical xtalk.  The output of the
01099   // routine is the number of xtalk pes generated in the xtalk pixel.
01100   //
01101   // For the configuration described by:
01102   // hp = hit pixel
01103   // hs = hit spot
01104   // xp = xtalk pixel
01105   // the database is contacted to obtain the parameter:
01106   // k = <xtalk_charge/hit_charge>. 
01107   //
01108   // These k values are taken from pmt test stand data and have 
01109   // been tuned to CalDet beam data, via the use of the global 
01110   // scaling parameters gOpticalXTScale, gOpticalXTOffset.
01111   //
01112   // Optical xtalk is stochastic in nature. Most of the time no charge is
01113   // produced. When xtalk is seen the charge corresponds to ~1pe.  
01114   // The average charge in the xtalk pixel is given by:
01115   //       average_xtalk_charge = k*hpe
01116   // with  k = <xtalk_charge/hit_charge>
01117   //
01118   // The number of PEs from xtalk is then generated by sampling a Poisson
01119   // distribution with a mean = average_xtalk_charge.
01120 
01121   // Can't work with nothing.
01122   if(hpe<=0.0) return 0.0;
01123   
01124   // Get a k value
01125   double k = RetrieveOpticalXtalkFraction(hp, hs, xp);
01126   // check k before going on
01127   if(k<=0.0) return 0.0; 
01128   
01129   // rescale k
01130   k=(k+gOpticalXTOffset)*gOpticalXTScale;
01131   
01132   // check k before going on
01133   if(k<=0.0) return 0.0; 
01134   
01135   // generate xtalk pes
01136   Int_t xpe = fRandom->Poisson(k*hpe);
01137   if(xpe<0) xpe=0; // some error from fRandom?
01138   
01139   return xpe;
01140   
01141 }

virtual Float_t SimPmtUTM16::GetDynodeGain (  )  const [inline, protected, virtual]

Reimplemented from SimPmt.

Definition at line 61 of file SimPmtUTM16.h.

00061 { return 0.8; };

virtual Pmt_t SimPmtUTM16::GetType (  )  const [inline, virtual]

Reimplemented from SimPmt.

Definition at line 56 of file SimPmtUTM16.h.

References SimPmt::kM16.

00056 { return SimPmt::kM16; };

bool SimPmtUTM16::InitSECValues (  )  [private]

Definition at line 770 of file SimPmtUTM16.cxx.

References PlexPixelSpotId::AsString(), fAverageSECValues, FractionalWidthToSEC(), fSECValues, SimPmt::fTube, SimPmt::GetGainAndWidth(), gkDefaultSECValues, Msg::kDebug, Msg::kError, Msg::kVerbose, MAXMSG, and MSG.

Referenced by Reset().

00771 {
00772   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00773   //
00774   // This routine fills the arrays fSECValues, fAverageSECValues
00775   // It should be called once when the object is made and again
00776   // whenever the gain and width information change for any pixel&spot
00777   // on this tube.
00778   
00779   MSG("DetSim",Msg::kDebug) << "InitSECValues() Tube " << fTube.AsString() << endl;
00780 
00781   static const double default_gain_last=gkDefaultSECValues[3]
00782     *gkDefaultSECValues[4]
00783     *gkDefaultSECValues[5]
00784     *gkDefaultSECValues[6]
00785     *gkDefaultSECValues[7]
00786     *gkDefaultSECValues[8]
00787     *gkDefaultSECValues[9]
00788     *gkDefaultSECValues[10]
00789     *gkDefaultSECValues[11];
00790 
00791   // rescale secs based on gains and widths taken from db
00792   // loop over pixels
00793   for(int pxit=0; pxit<16; pxit++){
00794     const int pixel=pxit+1;
00795     // loop over fiber spots on a pixel
00796     double ave_width=0.0;
00797     double ave_gain=0.0;
00798     for(int spit=0; spit<8; spit++){
00799       const int spot=spit+1;
00800       double pixel_gain, fractional_width;
00801       GetGainAndWidth(pixel,spot,pixel_gain,fractional_width);
00802 
00803       const double default_fractional_width=0.5;
00804 
00805       if(fractional_width<0.05 || fractional_width>5.0){
00806         // a crazy value for the spe width
00807         MAXMSG("DetSim", Msg::kError,50)
00808           <<"InitSECValues: GetGainAndWidth("<<pixel<<", "<<spot<<")"
00809           <<" returned fractional width " << fractional_width
00810           <<" which is outside the reasonable range!\n"
00811           <<" Will use the width "<< default_fractional_width
00812           <<" in further computations."<<endl;
00813         fractional_width=default_fractional_width;
00814       }
00815       const double sec = FractionalWidthToSEC(fractional_width);
00816       double width = fractional_width * pixel_gain;
00817       ave_width+=(width*width);
00818       ave_gain+=pixel_gain;
00819     // The following code sets the emmission coefficient of
00820     // the first 3 stages to sec and then rescales the 
00821     // default coefficients for the last 9 stages
00822     // so that the gain of the tube is equal to pixel_gain.
00823     // Since the width is largely determined by the emission
00824     // coefficients of the first 3 stages this procedure will result
00825     // in a simulation with the requested width and mean for the 1pe peak.
00826 
00827     // first 3 stages determine fractional width
00828       fSECValues[pxit][spit][0]=sec;
00829       fSECValues[pxit][spit][1]=sec;
00830       fSECValues[pxit][spit][2]=sec;
00831 
00832       const double pixel_gain_first=sec*sec*sec;
00833       // for other stages
00834       const double pixel_gain_last = pixel_gain/pixel_gain_first;
00835       const double rescale_factor = 
00836         pow((pixel_gain_last/default_gain_last),1.0/9.0);
00837       for(int stage=3; stage<12; stage++) {
00838         fSECValues[pxit][spit][stage] = 
00839           rescale_factor*gkDefaultSECValues[stage];
00840       }
00841     } // done with loop over spots
00842 
00843     // calculate average sec values
00844 
00845     // Now, ave_gain holds the sum of the gains for the 8 spots.
00846     // I divide by 8 to get the average.
00847     ave_gain/=8.0;
00848 
00849     // Now, ave_width holds the sum of the (fractional_width)^2 for
00850     // each spot. I then divide by 8 to average and take the square root.
00851     ave_width/=8.0;
00852     ave_width=sqrt(ave_width);
00853 
00854     // Then, the average secondary emission coeffiecient is computed from
00855     // ave_width, but needs to be converted back to a fractional width:
00856     double ave_sec=FractionalWidthToSEC(ave_width/ave_gain);
00857 
00858     // The first sec values of the first 3 stages are set to ave_sec.
00859     fAverageSECValues[pxit][0]=ave_sec;
00860     fAverageSECValues[pxit][1]=ave_sec;
00861     fAverageSECValues[pxit][2]=ave_sec;
00862 
00863     // The sec values of the last 9 stages are scaled to as to produce
00864     // an overall gain of ave_gain.
00865     const double ave_pixel_gain_first=ave_sec*ave_sec*ave_sec;
00866     const double ave_pixel_gain_last=ave_gain/ave_pixel_gain_first;
00867     // Figure out a rescaling factor.  The power of 1/9
00868     // is needed since we will end up taking this number to the 9th
00869     // power when computing the gain and/or simulating the dynode 
00870     // cascade.
00871     const double ave_rescale_factor = 
00872       pow((ave_pixel_gain_last/default_gain_last),1.0/9.0);
00873     for(int stage=3; stage<12; stage++) {
00874       fAverageSECValues[pxit][stage] = 
00875         ave_rescale_factor*gkDefaultSECValues[stage];
00876     }
00877     // done with rescaling business for this pixel
00878   }// done with loop over pixels
00879 
00880   // print out SEC values
00881   MSG("DetSim", Msg::kVerbose)<<"Printing SEC values: "<<endl;
00882   for(int pxit=0; pxit<16; pxit++){
00883     for(int spit=0; spit<8; spit++){
00884       MSG("DetSim", Msg::kVerbose)
00885         <<"Stage SECs for pixel "<<pxit+1<<", and spot "<<spit+1<<endl;
00886       for(int stit=0; stit<12; stit++){
00887       MSG("DetSim", Msg::kVerbose)
00888         <<"stage "<<stit+1<<" : "<<fSECValues[pxit][spit][stit]<<"   : (avg="
00889         <<fAverageSECValues[pxit][stit]<<")"<<endl;
00890       }
00891     }
00892   }
00893 
00894   return true;
00895 }

double SimPmtUTM16::LastStageChargeNonLinearity ( double  q,
double  gain,
double  lastsec 
) [static, private]

Definition at line 981 of file SimPmtUTM16.cxx.

References M16ChargeNonLinearity().

Referenced by OneHitSimDynodeChain().

00983 {
00984 
00985   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00986   //
00987   // This routine models pmt non linearity as being due to a change in
00988   // the secondary emission of the last dynode.  The data taken from
00989   // the M16 teststands and CalDet LI shows the amount of charge 
00990   // measured vs. the number of input PEs over the range of 0-300 input pes.
00991   //
00992   // The light level is determined by measuring the light level in the linear 
00993   // region of the tube (in practice ~30pes ), assuming linearity there, 
00994   // and then scaling by the known transmission values of the 
00995   // neutral density filters used in the stand or the PIN measurments.
00996   //
00997   // Qualitatively, for light levels > 50 pes a decrease in the amount of
00998   // charge per input pe is seen.  This decrease is treated as being due
00999   // to a change in the gain of the pmt.  
01000   
01001   // This routine uses the parameterization above to model the change
01002   // in gain by a deviation of the secondary emission coefficient of the last
01003   // dynode stage from its nominal value.  The parametrization is defined
01004   // in the member function M16ChargeNonLinearity.
01005 
01006   // The input parameters are:
01007   // q = charge input to the last stage in units of # electrons
01008   // gain = nominal gain of this pmt (unitless)
01009   // lsec = nominal sec of the last stage (unitless)
01010   
01011   double result=lsec;
01012   const double provisional_charge=q*lsec;
01013   const double frac_gain_change = M16ChargeNonLinearity(provisional_charge);
01014   result = (1.0+frac_gain_change)*lsec;
01015   
01016   return result;
01017 }

double SimPmtUTM16::LastStageNonLinearity ( double  q,
double  gain,
double  lastsec 
) [static, private]

Definition at line 897 of file SimPmtUTM16.cxx.

References Msg::kError, M16NonLinearity(), and MSG.

00898 {
00899   MSG("DetSim",Msg::kError) << "YOU ARE CALLING A DEPRECATED FUNCTION!"<<endl;
00900 
00901   //*** DEPRECATED FUNCTION ****   Mike Kordosky -- April 25, 2005
00902 
00903   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00904   //
00905   // This routine models pmt non linearity as being due to a change in
00906   // the secondary emission of the last dynode.  The data taken from
00907   // the M16 teststands shows the amount of charge measured vs. the number
00908   // of input PEs over the range of 0-300 input pes.  The light level
00909   // is determined by measuring the light level in the linear 
00910   // region of the tube (in practice ~30pes ) and then scaling by the 
00911   // known transmission values of the neutral density filters used 
00912   // in the stand.
00913   //
00914   // Qualitatively, for light levels > 50 pes a decrease in the amount of
00915   // charge per input pe is seen.  This decrease is treated as being due
00916   // to a change in the gain of the pmt.  The data is then processed to 
00917   // show the fractional change in gain vs. the number of input pes.
00918   // That data is then parameterized as:
00919   // (a) Flat for pe<=50
00920   // (b) 3rd order polynomial for 50<pe<=300
00921   // (c) For pe>300 the value at pe=300 is used
00922 
00923   // This routine uses the parameterization above to model the change
00924   // in gain by a deviation of the secondary emission coefficient of the last
00925   // dynode stage from its nominal value.  The parametrization is defined
00926   // in the member function M16NonLinearity.
00927 
00928   // The input parameters are:
00929   // q = charge input to the last stage
00930   // gain = nominal gain of this pmt
00931   // lsec = nominal sec of the last stage
00932   
00933   double result=lsec;
00934   const double approx_npe=q/(gain/lsec);
00935   if(approx_npe>50.0 && approx_npe<=300.0){
00936     const double frac_gain_change=M16NonLinearity(approx_npe);
00937     result = (1.0+frac_gain_change)*lsec;
00938   }
00939   else if(approx_npe>300.0){
00940     const double frac_gain_change=M16NonLinearity(300.0);
00941     result = (1.0+frac_gain_change)*lsec;
00942   }
00943 
00944   return result;
00945 }

double SimPmtUTM16::M16ChargeNonLinearity ( double  pe  )  [static, private]

Definition at line 1019 of file SimPmtUTM16.cxx.

References Munits::e_SI, SimPmt::fsVaGain, Msg::kDebug, and MSG.

Referenced by LastStageChargeNonLinearity().

01020 { 
01021   // Author: Mike Kordosky : kordosky@hep.utexas.edu
01022   //
01023   // Returns the fractional change in gain
01024   // due to pmt nonlinearity.
01025   //
01026   // Input: the charge, in # of electron units, at final stage
01027   //        typically this is a provisional charge found from the 
01028   //        result at the second to last stage multipled by the nominal
01029   //        gain of the last stage
01030   //
01031   // Parameterization taken from CalDet 2003 LI
01032   // thanks Anatael Cabrera
01033   //
01034   // See SimPmtUTM16::LastStageChargeNonLinearity
01035   // for more details.
01036   //
01037   
01038   // default fsVaGain = 0.375/Munits::fC
01039   // Munits::e_SI = 1.602176462e-19 C
01040   const double adc_per_electron = fsVaGain*Munits::e_SI;
01041 
01042   // calculate provisional adcs
01043   // need this because parameterization given in terms of VA ADCs
01044   // but we haven't digitized yet!
01045   const double provisional_adc = nelectrons*adc_per_electron;
01046   
01047   // calculate gain dip
01048   // gain_dip: change of gain at last stage owing to non linearity
01049   //           negative # means a loss in tube gain
01050   //           the gain of the last stage should be modeled as
01051   //           gain = nominal*(1+gain_dip)
01052   //
01053   static const double p[2]={0.0166, -4.51e-6};
01054   
01055   double gain_dip = p[0] + p[1]*provisional_adc;
01056 
01057   const double linear_adc = 3600;//value below which PMTs are linear
01058   const double high_adc = 60*300;//upper limit of the parameterization ~300 PEs
01059   if(provisional_adc<=linear_adc) { 
01060     MSG("DetSim",Msg::kDebug) << " Tube in linear region: "
01061                               <<" provisional adc = "<<provisional_adc
01062                               <<"  < linear_adc = "<<linear_adc<<"\n"
01063                               <<" gain_dip set to zero "<<endl;
01064     gain_dip=0.0;
01065   }
01066   else if(provisional_adc>high_adc){
01067     gain_dip = p[0] + high_adc*p[1];
01068     MSG("DetSim",Msg::kDebug) << " Tube outside of parameterized region: "
01069                               <<" provisional adc = "<<provisional_adc
01070                               <<"  < high_adc = "<<high_adc<<"\n"
01071                               <<" gain_dip set to value at"
01072                               <<high_adc<<" ADCs : "<<gain_dip<<endl;
01073     
01074   }
01075 
01076   MSG("DetSim",Msg::kDebug) << "[nelectrons, provisional adc, %gain_dip] : [ " 
01077                             << nelectrons<<" , "
01078                             << provisional_adc<<" , "
01079                             << gain_dip*100.0 <<" ]"<<endl;
01080 
01081 
01082   return gain_dip; 
01083 
01084 }

double SimPmtUTM16::M16NonLinearity ( double  pe  )  [static, private]

Definition at line 947 of file SimPmtUTM16.cxx.

References Msg::kError, and MSG.

Referenced by LastStageNonLinearity().

00948 { 
00949   MSG("DetSim",Msg::kError) << "YOU ARE CALLING A DEPRECATED FUNCTION!"<<endl;
00950   
00951   //***DEPRECATED FUNCTION***  Mike Kordosky -- April 25, 2005
00952   //
00953   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00954   //
00955   // Returns the fractional change in gain
00956   // due to pmt nonlinearity.
00957   // Parameterization taken from m16 teststand data.
00958   // See SimPmtUTM16::LastStageNonLinearity
00959   // for more details.
00960   //
00961   // This function is to be applied for 50<pe<300.
00962   
00963   //  static const double p[4]={5.0E-2,
00964   //                        -1.4E-3,
00965   //                        5.6E-6,
00966   //                        7.5E-9};
00967 
00968   // M. Kordosky - March 22, 2005
00969   // sign error in 3rd degree term... I thought I fixed it long ago!!
00970   // thanks to Anatael for illuminating this problem
00971   //
00972   static const double p[4]={5.0E-2,
00973                             -1.4E-3,
00974                             5.6E-6,
00975                             -7.5E-9};
00976   
00977   return p[0] + p[1]*pe + p[2]*pe*pe + p[3]*pe*pe*pe;
00978 }

bool SimPmtUTM16::OneHitSimDynodeChain ( Int_t  hp,
Int_t  hs,
Float_t  hpe,
std::vector< double > &  qvec,
std::vector< double > &  tvec 
) [private]

Definition at line 449 of file SimPmtUTM16.cxx.

References Munits::e_SI, fAverageSECValues, SimPmt::fRandom, fSECValues, SimPmt::fsPmtDoChargeCrosstalk, SimPmt::fsPmtDoNonlinearity, SimPmt::fsPmtScaleChargeCrosstalk, gElectricXTOffset, gElectricXTScale, gTransitTimeJitterNS, gTransitTimeNS, Msg::kVerbose, LastStageChargeNonLinearity(), MSG, Munits::ns, PrintCharge(), RetrieveElectricXtalkFraction(), and ScaleDiagonalAdjacentChargeCrosstalk().

Referenced by SimulateCharges().

00452 {
00453   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00454   //
00455   // Starting with pes on one pixel, this routine
00456   // simulates the secondary emission process at each dynode
00457   // stage.  At each stage after the first, prior to the calculation
00458   // of secondary electrons, this routine generates xtalk by leaking
00459   // some fraction of the charge in the hit pixel (at the current stage)
00460   // into the other pixels of the tube.  The leakage is done in a stochasitic
00461   // fashion with the average fraction of charge leaked taken from a database
00462   // table. After the simulation of charge leakage the number of secondary
00463   // electrons for each pixel is calculated by sampling a poisson distribution
00464   // with a mean given by:
00465   // (charge at current stage and pixel)
00466   // *(secondary emission coefficient at current stage and pixel)
00467   //
00468   // The charges created in this procedure are ADDED to what is already
00469   // in qvec.
00470 
00471   // NOTE:  Much of the rescaling work can probably be moved outside 
00472   // of this function so as to limit the number of times it must be done.
00473   // MAK 2003.05.12 -- now done
00474 
00475   // If the function was passed zero pes then quit.
00476   if(hpe<=0.0) return true;
00477 
00478   // Set some array indices
00479   const int hpit=hp-1;
00480   const int hsit=hs-1;
00481   MSG("DetSim", Msg::kVerbose)<<"Calculating gains."<<endl;
00482   // calculate the gains at each pixel
00483   double gains[16];
00484   for(int pxit=0; pxit<16; pxit++){
00485     gains[pxit]=1.0;
00486     for(int stage=0; stage<12; stage++){
00487       // for the hit pixel, we can use the exact gain 
00488       // since the spot is known
00489       if(pxit==hpit) gains[pxit]*=fSECValues[hpit][hsit][stage];
00490       // for xtalk pixels use an average gain
00491       else gains[pxit]*=fAverageSECValues[pxit][stage];
00492     }
00493   }
00494 
00495   // Now simulate dynode amplification process
00496   // This includes electric xtalk and pmt nonlinearity
00497   
00498   // an array to hold the charge on each pixel
00499   // gets updated as this routine proceeds
00500   unsigned int qarray[16]={0};
00501   
00502   // look up electric xtalk fractions
00503   // hold the fractions in xtfrac 
00504   MSG("DetSim", Msg::kVerbose)<<"Looking up xtalk fractions."<<endl;
00505   double xtfrac[16]={0.0};
00506   for(int pxit=0; pxit<16; pxit++){
00507     if(pxit==hpit) continue;
00508     const double got_xtfrac=RetrieveElectricXtalkFraction(hp,hs,pxit+1);
00509     if(got_xtfrac<0.0) xtfrac[pxit]=0.0;
00510     else xtfrac[pxit]=got_xtfrac;
00511   }
00512 
00513   // first stage: amplify the hit pes
00514   const unsigned int begin_charge=
00515     (unsigned int) ( fRandom->PoissonD(hpe*fSECValues[hpit][hsit][0]) );
00516   qarray[hpit]=begin_charge;
00517   MSG("DetSim", Msg::kVerbose)<<"Charge after stage 1"<<endl;
00518   PrintCharge(qarray);
00519   // loop through the last 11 stages
00520   // at each stage, leak some electrons for electric xtalk
00521   // then do dynode multiplication
00522   // NOTE: leakage prior to the first dynode looks like
00523   // what we normally call optical xtalk... that is why this
00524   // algorithm starts on the second stage.
00525   for(int stage=1; stage<12; stage++){ 
00526     // first leak charge for xtalk
00527     if(fsPmtDoChargeCrosstalk){
00528       for(int pxit=0; pxit<16; pxit++){
00529         if(pxit==hpit) continue; // pixels do not xtalk to themselves
00530         if(qarray[hpit]<=0.0) break; // can't work with nothing
00531         // create some xtalk charge
00532         double working_xt_fraction
00533           =(xtfrac[pxit]+gElectricXTOffset)*gElectricXTScale/11.0;
00534         // Scale it by user.
00535         working_xt_fraction *= fsPmtScaleChargeCrosstalk;
00536         
00537   //Below code added by G. Pawloski to scale diagonal and adjacent Xtalk separately 
00538         working_xt_fraction *= ScaleDiagonalAdjacentChargeCrosstalk(pxit,hpit);
00539 
00540         qarray[pxit]+=(unsigned int)(fRandom->PoissonD(working_xt_fraction*qarray[hpit]));
00541       }
00542     }
00543     // amplify charges at this stage, including any charge from xtalk
00544     for(int pxit=0; pxit<16; pxit++){
00545       if(qarray[pxit]<=0) continue;// can't work with nothing
00546       double current_sec=0.0;
00547       if(pxit==hpit) current_sec=fSECValues[hpit][hsit][stage];
00548       else current_sec=fAverageSECValues[hpit][stage];
00549       // correct for non linearity on the last stage
00550       if(stage==11 && fsPmtDoNonlinearity){
00551         current_sec=
00552           LastStageChargeNonLinearity(qarray[pxit],gains[pxit],current_sec);
00553       }
00554       // simulate charge amplification at this stage
00555       qarray[pxit]=(unsigned int)(fRandom->PoissonD(qarray[pxit]*current_sec));
00556 
00557     }
00558     MSG("DetSim", Msg::kVerbose)<<"Charge after stage "<<(stage+1)<<endl;
00559     PrintCharge(qarray);
00560   }
00561 
00562 
00563   // dynode/anode charge fraction
00564   // Teststand data shows dynode/anode ~ 35%-50%
00565   const double dynode_fraction=0.425;
00566   // ADD created charges into ouput
00567   for(int pxit=0; pxit<16; pxit++){
00568     // sum charges for dynode output
00569     qvec[16]+=(qarray[pxit]/dynode_fraction)*Munits::e_SI;
00570     // convert charge from # of electrons to C
00571     qvec[pxit]+=qarray[pxit]*Munits::e_SI; // add charge in C
00572   }
00573 
00574   // The timing:
00575   // Hamamatsu says transit time 8.8 ns with a jitter of 180 ps
00576   // these numbers are for single pes.
00577   // The input charge drives the generation of charge in the tube, 
00578   // so simulate for it and apply the result to all pixels.
00579 
00580   // what is the common currency for times in DetSim?
00581   // I'm guessing seconds based on the above comment regarding C vs. fC.
00582   const double transit_time = 
00583     fRandom->Gaus(gTransitTimeNS, 
00584                   (gTransitTimeJitterNS)/sqrt(hpe))*Munits::ns;
00585   for(int pxit=0; pxit<17; pxit++){
00586     tvec[pxit]+=transit_time;
00587   }
00588 
00589   return true;
00590 }

void SimPmtUTM16::Print ( Option_t *  option = ""  )  const [virtual]

Reimplemented from SimPmt.

Definition at line 226 of file SimPmtUTM16.cxx.

References PlexPixelSpotId::AsString(), Munits::fC, SimPmt::fTube, SimPmt::GetCharge(), SimPmt::GetPe(), SimPmt::GetPeXtalk(), and it.

00227 {
00228   printf("  SimPmtUTM16::Print()  -- Tube: %s\n",fTube.AsString("t"));
00229   
00230   std::string opt = option;
00231   Bool_t spots = (opt.find("spot") != std::string::npos);;
00232 
00233   for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
00234     int ibucket = it.BucketId();
00235     
00236   
00237     printf("  Time bucket %d:\n",ibucket);
00238 
00239     // output:
00240     // pixels:  16 15 14 13   fibres:  8 7 6
00241     //          12 11 10  9             5 4
00242     //           8  7  6  5            3 2 1
00243     //           4  3  2  1
00244   
00245     if(spots) { 
00246       printf("  Photoelectrons(Npe)\n");
00247       for(int row=0; row < 4; row++) {
00248         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",
00249                GetPe(16-row*4,8,ibucket),  GetPe(16-row*4,7,ibucket),  GetPe(16-row*4,6,ibucket),
00250                GetPe(15-row*4,8,ibucket),  GetPe(15-row*4,7,ibucket),  GetPe(15-row*4,6,ibucket),
00251                GetPe(14-row*4,8,ibucket),  GetPe(14-row*4,7,ibucket),  GetPe(14-row*4,6,ibucket),
00252                GetPe(13-row*4,8,ibucket),  GetPe(13-row*4,7,ibucket),  GetPe(13-row*4,6,ibucket));
00253 
00254         printf("     %5.0f %5.0f        %5.0f %5.0f        %5.0f %5.0f        %5.0f %5.0f\n",
00255                GetPe(16-row*4,5,ibucket),  GetPe(16-row*4,4,ibucket),
00256                GetPe(15-row*4,5,ibucket),  GetPe(15-row*4,4,ibucket),
00257                GetPe(14-row*4,5,ibucket),  GetPe(14-row*4,4,ibucket),
00258                GetPe(13-row*4,5,ibucket),  GetPe(13-row*4,4,ibucket));
00259     
00260         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",
00261                GetPe(16-row*4,3,ibucket),  GetPe(16-row*4,2,ibucket),  GetPe(16-row*4,1,ibucket),
00262                GetPe(15-row*4,3,ibucket),  GetPe(15-row*4,2,ibucket),  GetPe(15-row*4,1,ibucket),
00263                GetPe(14-row*4,3,ibucket),  GetPe(14-row*4,2,ibucket),  GetPe(14-row*4,1,ibucket),
00264                GetPe(13-row*4,3,ibucket),  GetPe(13-row*4,2,ibucket),  GetPe(13-row*4,1,ibucket));
00265            
00266         printf("\n");
00267       }
00268     }
00269   
00270   
00271     printf("\n        PEs                   PE with Xtalk         Charges (fC)\n");
00272     for(int row=0; row < 4; row++) {
00273       printf("    %4d %4d %4d %4d",
00274              (int)GetPe(16-row*4,0,ibucket), (int)GetPe(15-row*4,0,ibucket),
00275              (int)GetPe(14-row*4,0,ibucket), (int)GetPe(13-row*4,0,ibucket));
00276     
00277       printf("    %4d %4d %4d %4d",
00278              (int)GetPeXtalk(16-row*4,0,ibucket), (int)GetPeXtalk(15-row*4,0,ibucket),
00279              (int)GetPeXtalk(14-row*4,0,ibucket), (int)GetPeXtalk(13-row*4,0,ibucket));
00280 
00281       printf("    %6.0f %6.0f %6.0f %6.0f\n",
00282              GetCharge(16-row*4,ibucket)/Munits::fC, GetCharge(15-row*4,ibucket)/Munits::fC,
00283              GetCharge(14-row*4,ibucket)/Munits::fC, GetCharge(13-row*4,ibucket)/Munits::fC);
00284     }
00285   
00286   }
00287   
00288   printf("\n");
00289 }

void SimPmtUTM16::PrintCharge ( unsigned int *  q  )  [private]

Definition at line 1157 of file SimPmtUTM16.cxx.

References Msg::kVerbose, and MSG.

Referenced by OneHitSimDynodeChain().

01157                                                  {
01158   static const MsgFormat ifmt("%10d");
01159   
01160   MSG("DetSim", Msg::kVerbose)
01161     <<"================================================================================\n"
01162     <<ifmt(qarray[3])<<" "<<ifmt(qarray[2])<<" "<<ifmt(qarray[1])<<" "<<ifmt(qarray[0])
01163     <<"\n"
01164     <<ifmt(qarray[7])<<" "<<ifmt(qarray[6])<<" "<<ifmt(qarray[5])<<" "<<ifmt(qarray[4])
01165     <<"\n"
01166     <<ifmt(qarray[11])<<" "<<ifmt(qarray[10])<<" "<<ifmt(qarray[9])<<" "<<ifmt(qarray[8])
01167     <<"\n"
01168     <<ifmt(qarray[15])<<" "<<ifmt(qarray[14])<<" "<<ifmt(qarray[13])<<" "<<ifmt(qarray[12])
01169     <<"\n"
01170     <<"================================================================================\n"
01171 
01172     <<endl;
01173 }

void SimPmtUTM16::Reset ( const VldContext context  )  [virtual]

Reimplemented from SimPmt.

Definition at line 185 of file SimPmtUTM16.cxx.

References fSECValuesBuilt, SimPmt::fsRebuildGainMap, gCrosstalkTable, InitSECValues(), Msg::kFatal, MSG, and SimPmtM16CrosstalkTable::Reset().

00186 {
00187   // Author: N. Tagg
00188   //
00189   // Resets the PMT for a new event at a new context.
00190   // Ensures the correct DB is used.
00191 
00192   // Reset the PMT data.
00193   SimPmt::Reset(context);
00194 
00195   // Pulls new xtalk data if available, or 
00196   // gets it for the first time if it doesn't yet exist.
00197   if(gCrosstalkTable) {
00198     gCrosstalkTable->Reset(context);
00199   } else {
00200     gCrosstalkTable = new SimPmtM16CrosstalkTable(context);
00201   }
00202 
00203   // R1.17 Reboot the model:
00204   if(fsRebuildGainMap || (!fSECValuesBuilt)) 
00205     fSECValuesBuilt = InitSECValues();
00206 
00207   if(!fSECValuesBuilt){
00208     MSG("DetSim", Msg::kFatal)<<"Could not initialize SEC values!"<<endl;
00209     assert(0); // make program bomb out... is this the right thing to do?    
00210   }
00211 }

double SimPmtUTM16::RetrieveElectricXtalkFraction ( Int_t  hp,
Int_t  hs,
Int_t  xp 
) const [private]

Definition at line 681 of file SimPmtUTM16.cxx.

References gCrosstalkTable, SimPmtM16Crosstalk::GetElecFrac(), SimPmtM16CrosstalkTable::GetRow(), Msg::kError, Msg::kFatal, Msg::kVerbose, and MSG.

Referenced by OneHitSimDynodeChain().

00684 {
00685   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00686   //
00687   // This routine should talk to the database and get the appropriate constant.
00688   // It can employ this object's fContext member and 
00689   // the arguments to this function.
00690   
00691   double result=0.0;
00692   if(gCrosstalkTable){
00693     const SimPmtM16Crosstalk* row = gCrosstalkTable->GetRow(hp,hs,xp);
00694     if(row){
00695       result = row->GetElecFrac();
00696       static const MsgFormat ffmt("%9.3e");
00697       MSG("DetSim", Msg::kVerbose)
00698         <<"RetrieveElectricXtalkFraction("<<hp<<", "<<hs<<", "<<xp<<") "
00699         <<"will return "<<ffmt(result)<<" ."<<endl;
00700     }
00701     else{
00702       MSG("DetSim", Msg::kError)
00703         <<"RetrieveElectricXtalkFraction "
00704         <<"could not find a db row for:\n"
00705         <<"hit_pixel: "<< hp
00706         <<", hit_spot: "<<hs
00707         <<", xtalk_pixel: "<<xp
00708         <<"."<<endl;
00709     }
00710   }
00711   else{
00712     MSG("DetSim", Msg::kFatal)
00713       <<"Somehow RetrieveElectricXtalkFraction was called with "
00714       "gCrosstalkTable==0!\nThis is a fatal error.";
00715       assert(0); // bomb out
00716   }
00717 
00718   return result; // a return value <= 0.0 will result in no xtalk
00719 }

double SimPmtUTM16::RetrieveOpticalXtalkFraction ( Int_t  hp,
Int_t  hs,
Int_t  xp 
) const [private]

Definition at line 721 of file SimPmtUTM16.cxx.

References SimPmt::fsPmtScaleOpticalCrosstalk, gCrosstalkTable, SimPmtM16Crosstalk::GetOptFrac(), SimPmtM16CrosstalkTable::GetRow(), Msg::kError, Msg::kFatal, Msg::kVerbose, MSG, and ScaleDiagonalAdjacentOpticalCrosstalk().

Referenced by GenOpticalCrosstalk().

00724 {
00725   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00726   //
00727   // This routine should talk to the database and get the appropriate constant.
00728   // It can employ this object's fContext member and 
00729   // the arguments to this function.
00730 
00731   double result=0.0;
00732   if(gCrosstalkTable){
00733     const SimPmtM16Crosstalk* row = gCrosstalkTable->GetRow(hp,hs,xp);
00734     if(row){
00735       result = row->GetOptFrac();
00736       static const MsgFormat ffmt("%9.3e");
00737       MSG("DetSim", Msg::kVerbose)
00738         <<"RetrieveOpticalXtalkFraction("<<hp<<", "<<hs<<", "<<xp<<") "
00739         <<"will return "<<ffmt(result)<<" ."<<endl;
00740     }
00741     else{
00742       MSG("DetSim", Msg::kError)
00743         <<"RetrieveOpticalXtalkFraction "
00744         <<"could not find a db row for:\n"
00745         <<"hit_pixel: "<< hp
00746         <<", hit_spot: "<<hs
00747         <<", xtalk_pixel: "<<xp
00748         <<"."<<endl;
00749     }
00750   }
00751   else{
00752     MSG("DetSim", Msg::kFatal)
00753       <<"Somehow RetrieveOpticalXtalkFraction was called with "
00754       "gCrosstalkTable==0!\nThis is a fatal error.";
00755       assert(0); // bomb out
00756   }
00757 
00758 
00759   //Below code added by G. Pawloski to scale diagonal and adjacent Xtalk separately 
00760 //  return result * fsPmtScaleOpticalCrosstalk; // a return value <= 0.0 will result in no xtalk
00761 
00762   result *= fsPmtScaleOpticalCrosstalk;
00763 
00764         result *= ScaleDiagonalAdjacentOpticalCrosstalk(xp-1,hp-1);
00765   return result; // a return value <= 0.0 will result in no xtalk
00766 
00767 }

double SimPmtUTM16::ScaleDiagonalAdjacentChargeCrosstalk ( Int_t  pxit,
Int_t  hpit 
) const [private]

Definition at line 595 of file SimPmtUTM16.cxx.

References adjacentPixels, diagonalPixels, SimPmt::fsPmtScaleAdjacentChargeCrosstalk, and SimPmt::fsPmtScaleDiagonalChargeCrosstalk.

Referenced by OneHitSimDynodeChain().

00596 {
00597   // Author: Greg Pawloski : pawloski@fnal.gov
00598   //
00599   // This routine applies a separate scaling factor for diagonal and adjacent
00600   // electric xtalk.  Routine check if xtalk pixel (pxit) is diagonal and adjacent to
00601   // hit pixel (hpit)
00602   //
00603   //   Pixels:
00604   //   3  2  1  0
00605   //   7  6  5  4
00606   //   11 10 9  8
00607   //   15 14 13 12
00608   // There might be a smarter way of doing this but using a lookup table is less prone to logic errors
00609 
00610  double scale=1.0;
00611  if(
00612         pxit < 0
00613      || pxit > 15
00614      || hpit < 0
00615      || hpit > 15
00616    )
00617   {
00618    return(1.0);
00619   }
00620  bool isAdjacentPixels = (adjacentPixels[hpit].find(pxit) != adjacentPixels[hpit].end());
00621  bool isDiagonalPixels = (diagonalPixels[hpit].find(pxit) != diagonalPixels[hpit].end());
00622 
00623  if(isAdjacentPixels)
00624  {
00625    scale = fsPmtScaleAdjacentChargeCrosstalk;
00626  }
00627  else if(isDiagonalPixels)
00628  {
00629    scale = fsPmtScaleDiagonalChargeCrosstalk;
00630  }
00631 
00632  return(scale);
00633 
00634 }

double SimPmtUTM16::ScaleDiagonalAdjacentOpticalCrosstalk ( Int_t  pxit,
Int_t  hpit 
) const [private]

Definition at line 636 of file SimPmtUTM16.cxx.

References adjacentPixels, diagonalPixels, SimPmt::fsPmtScaleAdjacentOpticalCrosstalk, and SimPmt::fsPmtScaleDiagonalOpticalCrosstalk.

Referenced by RetrieveOpticalXtalkFraction().

00637 {
00638   // Author: Greg Pawloski : pawloski@fnal.gov
00639   //
00640   // This routine applies a separate scaling factor for diagonal and adjacent
00641   // electric xtalk.  Routine check if xtalk pixel (pxit) is diagonal and adjacent to
00642   // hit pixel (hpit)
00643   //
00644   //   Pixels:
00645   //   3  2  1  0
00646   //   7  6  5  4
00647   //   11 10 9  8
00648   //   15 14 13 12
00649   // There might be a smarter way of doing this but using a lookup table is less prone to logic errors
00650 
00651  double scale=1.0;
00652  if(
00653         pxit < 0
00654      || pxit > 15
00655      || hpit < 0
00656      || hpit > 15
00657    )
00658   {
00659    return(1.0);
00660   }
00661  bool isAdjacentPixels = (adjacentPixels[hpit].find(pxit) != adjacentPixels[hpit].end());
00662  bool isDiagonalPixels = (diagonalPixels[hpit].find(pxit) != diagonalPixels[hpit].end());
00663 
00664  if(isAdjacentPixels)
00665  {
00666    scale = fsPmtScaleAdjacentOpticalCrosstalk;
00667  }
00668  else if(isDiagonalPixels)
00669  {
00670    scale = fsPmtScaleDiagonalOpticalCrosstalk;
00671  }
00672 
00673  return(scale);
00674 
00675 }

void SimPmtUTM16::SimulateCharges (  )  [protected, virtual]

Reimplemented from SimPmt.

Definition at line 380 of file SimPmtUTM16.cxx.

References SimPixelTimeBucket::AddCharge(), SimPmtTimeBucket::AddTotalCharge(), SimPmt::fNPixels, SimPmt::fNSpots, SimPmt::fTotalCharge, SimPixelTimeBucket::GetPEXtalk(), SimPmtTimeBucket::GetPixelBucket(), it, DigiSignal::kCrosstalk, Msg::kError, Msg::kVerbose, MSG, OneHitSimDynodeChain(), and SimPixelTimeBucket::SetCharge().

Referenced by SimulatePmt().

00381 {
00382 
00383   MSG("DetSim", Msg::kVerbose)
00384     <<"Calling SimPmtUTM16::SimulateCharges()"<<endl;
00385   
00386   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00387   //
00388   // basic structure taken from SimPmt::SimulateCharges()
00389   // this proceedure simulates charge xtalk in parallel
00390   // with the dynode amplification process
00391   //
00392   fTotalCharge = 0;
00393   // Iterate over all time buckets.
00394   MSG("DetSim", Msg::kVerbose)<<"Iterating over time buckets now."<<endl;
00395 
00396   int cntr=1;
00397   for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
00398     SimPmtTimeBucket& pmttb = it.Bucket();
00399     MSG("DetSim", Msg::kVerbose)<<"Working on time bucket: "<<cntr<<endl;
00400     cntr++;
00401     // zero charge in all pixels of this time bucket
00402     // original code did the equivalent
00403     // not sure if it is really needed
00404     for(int pix=1; pix<=fNPixels; pix++) pmttb.GetPixelBucket(pix).SetCharge(0);
00405     // Iterate over hit pixels in this time bucket.
00406     MSG("DetSim", Msg::kVerbose)<<"Iterating over hit pixels."<<endl;
00407     for(Int_t hitpix = 1; hitpix <= fNPixels; hitpix++) {
00408       MSG("DetSim", Msg::kVerbose)<<"Working on pixel: "<<hitpix<<endl;
00409       bool oksim=false;
00410       for(Int_t spot = 1; spot <= fNSpots; spot++) {
00411         const double hitpe = pmttb.GetPixelBucket(hitpix).GetPEXtalk(spot);
00412         MSG("DetSim", Msg::kVerbose)<<"Working on spot: "<<spot
00413                                   <<" got hitpe= "<<hitpe<<endl;
00414         if(hitpe<=0.0) continue;
00415         // these collections will hold the results of the simulation
00416         std::vector<double> qvec(fNPixels+1, 0.0); // charge
00417         std::vector<double> tvec(fNPixels+1, 0.0); // transit time
00418         MSG("DetSim", Msg::kVerbose)<<"Calling OneHitSimDynodeChain("
00419                                   <<hitpix<<", "<<spot<<", "<<hitpe
00420                                   <<", qvec, tvec)"<<endl;      
00421         oksim = OneHitSimDynodeChain(hitpix, spot, hitpe, qvec, tvec);
00422         // add the results to the pixel time buckets
00423         if(oksim){
00424           for(int outpix=1; outpix<=fNPixels; outpix++){
00425             const double generated_charge = qvec[outpix-1];
00426             MSG("DetSim", Msg::kVerbose)<<"Generated charge="<<generated_charge
00427                                       <<" on pixel "<<outpix
00428                                       <<endl;
00429             pmttb.GetPixelBucket(outpix).AddCharge(generated_charge);
00430             // if this is not the hit pixel, then any charge is from xtalk
00431             if(generated_charge>0 && outpix!=hitpix){
00432               // flip that xtalk bit
00433               pmttb.GetPixelBucket(outpix).
00434                 SetTruthBit(DigiSignal::kCrosstalk);
00435             }
00436           } 
00437           pmttb.AddTotalCharge(qvec[fNPixels]);
00438           fTotalCharge+=qvec[fNPixels];
00439           // what do I do with the transit time information?
00440         }
00441         else{
00442           MSG("DetSim", Msg::kError)<<"Error in OneHitSimDynodeChain.\nDid not place any charge into the time bucket.\nWill try to continue."<<endl;
00443         }
00444       } // spots
00445     } // pixels
00446    } // buckets.
00447 }

void SimPmtUTM16::SimulatePmt ( void   )  [protected, virtual]

Reimplemented from SimPmt.

Definition at line 213 of file SimPmtUTM16.cxx.

References SimPmt::CopyPEtoPEXtalk(), SimPmt::fsPmtDoDarkNoise, SimPmt::fsPmtDoOpticalCrosstalk, SimulateCharges(), SimPmt::SimulateDarkNoise(), and SimPmt::SimulateOpticalXtalk().

00214 {
00215   // Author: Mike Kordosky : kordosky@hep.utexas.edu
00216 
00217   if(fsPmtDoOpticalCrosstalk) SimulateOpticalXtalk();
00218   else CopyPEtoPEXtalk();  // A null operation.
00219 
00220   if(fsPmtDoDarkNoise) SimulateDarkNoise();
00221   SimulateCharges(); // internally uses fElecticXtalk, fNonLinearity
00222   
00223 }


Member Data Documentation

std::set<int> SimPmtUTM16::adjacentPixels[16] [private]
std::set<int> SimPmtUTM16::diagonalPixels[16] [private]
double SimPmtUTM16::fAverageSECValues[16][12] [private]

Definition at line 113 of file SimPmtUTM16.h.

Referenced by InitSECValues(), and OneHitSimDynodeChain().

Definition at line 72 of file SimPmtUTM16.h.

double SimPmtUTM16::fSECValues[16][8][12] [private]

Definition at line 110 of file SimPmtUTM16.h.

Referenced by InitSECValues(), and OneHitSimDynodeChain().

Definition at line 107 of file SimPmtUTM16.h.

Referenced by Reset().

double SimPmtUTM16::gElectricXTOffset = -0.001 [static, private]

Definition at line 101 of file SimPmtUTM16.h.

Referenced by OneHitSimDynodeChain().

double SimPmtUTM16::gElectricXTScale = 0.75 [static, private]

Definition at line 100 of file SimPmtUTM16.h.

Referenced by OneHitSimDynodeChain().

const double SimPmtUTM16::gkDefaultGain [static, private]
const double SimPmtUTM16::gkDefaultSECValues [static, private]
Initial value:
 {4.9407, 4.9407, 4.9407,
                                                    2.53997, 2.53997, 2.53997,
                                                    2.53997, 2.53997, 2.53997,
                                                    2.53997, 2.91747, 4.9407}

Definition at line 96 of file SimPmtUTM16.h.

Referenced by InitSECValues().

double SimPmtUTM16::gOpticalXTOffset = 0.0 [static, private]

Definition at line 99 of file SimPmtUTM16.h.

Referenced by GenOpticalCrosstalk().

double SimPmtUTM16::gOpticalXTScale = 1.75 [static, private]

Definition at line 98 of file SimPmtUTM16.h.

Referenced by GenOpticalCrosstalk().

double SimPmtUTM16::gTransitTimeJitterNS = 0.180 [static, private]

Definition at line 103 of file SimPmtUTM16.h.

Referenced by OneHitSimDynodeChain().

double SimPmtUTM16::gTransitTimeNS = 8.8 [static, private]

Definition at line 102 of file SimPmtUTM16.h.

Referenced by OneHitSimDynodeChain().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1