SimPmtM64.cxx

Go to the documentation of this file.

#include "SimPmtM64.h"
#include "SimPmtM64CrosstalkTable.h"
#include "SimPmtMaker.h"
#include "MessageService/MsgService.h"
#include "TMath.h"

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

SimPmtMakerProxy<SimPmtM64> gSimPmtM64Proxy("SimPmtM64");

ClassImp(SimPmtM64)

// Crosstalk database table:
SimPmtM64CrosstalkTable* SimPmtM64::fsCrosstalkTable = NULL;

SimPmtM64::SimPmtM64(PlexPixelSpotId tube, 
            VldContext context,
            TRandom* random ) :
  SimPmt(tube,context,random,64,1),
  fQieClock(context)
{
  // Constructor.
  
  // Single global reference to the crosstalk table.
  // This assumes that the crosstalk table has no validity 
  // changes; may need to change this in the future.
  //
  // Note that this is a global memory leak, too. Really should 
  // think about fixing it.
  if(fsCrosstalkTable == NULL) {
    fsCrosstalkTable = new SimPmtM64CrosstalkTable(context);
  };


  // Build cache of gains/widths.
  Double_t g, w;
  for(Int_t pix = 1; pix <= fNPixels; pix++) {      
    GetGainAndWidth(pix, 1, g, w);
    // Protect against crazy widths:
    if(w>0.8) w = 0.8;    

    // The secondary emmission ratio is directly related 
    // to the width of the 1pe peak.
    // In fact, this is 1/(w^2) where w is the fractional
    // width of the 1pe peak.
    Float_t secemm = 1.0/(w*w);  if(w>0.8) w = 0.8;
    
    fPixelGains[pix] = g; // Gain
    fPixelSecEmm[pix] = secemm; // secondary emmission ratio.
  };
}

void     SimPmtM64::Reset( const VldContext& context )
{
  // Reset the PMT for a new event.
  // Reset the clock:
  fQieClock.Reset(context);

  // Make sure database is up-to-date
  if(fsCrosstalkTable) fsCrosstalkTable->Reset(context);
  else fsCrosstalkTable = new SimPmtM64CrosstalkTable(context);
  // Reset the data:
  SimPmt::Reset(context);

}

Int_t    SimPmtM64::TimeToBucket(Double_t time)
{
  // Use the QIE clock to assign time buckets.
  return fQieClock.GetBucketForTime(time);
}

Double_t    SimPmtM64::BucketToStartTime(Int_t bucket)
{
  // Use the QIE clock to assign time buckets.
  return fQieClock.GetTimeOfBucket(bucket);
}

Double_t    SimPmtM64::BucketToStopTime(Int_t bucket)
{
  // Use the QIE clock to assign time buckets.
  return fQieClock.GetTimeOfBucket(bucket+1);
}


Float_t  SimPmtM64::GetBucketDuration( Int_t /* bucket */ )
{
  return fQieClock.fNDTick;
}

void SimPmtM64::Print(Option_t* /* option */) const
{
  printf("  SimPmtM64::Print()  -- Tube: %s\n",fTube.AsString("t"));
  
  SimPmtBucketIterator it(*(SimPmt*)this);
  printf("      Photoelectrons(Npe)                        Charge (fC)   \n");
  for( ; !it.End(); it.Next() ) {
    int ibucket = it.BucketId();
    printf("  ---------------------------Time bucket %d----------------------\n",ibucket);
    for(int row=0; row < 8; row++) {
      printf("  %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f",
             GetPe(row*8+1,1,ibucket),
             GetPe(row*8+2,1,ibucket),
             GetPe(row*8+3,1,ibucket),
             GetPe(row*8+4,1,ibucket),
             GetPe(row*8+5,1,ibucket),
             GetPe(row*8+6,1,ibucket),
             GetPe(row*8+7,1,ibucket),
             GetPe(row*8+8,1,ibucket));
      
      printf(" |  %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f",
             GetCharge(row*8+1,ibucket)/Munits::fC, 
             GetCharge(row*8+2,ibucket)/Munits::fC,
             GetCharge(row*8+3,ibucket)/Munits::fC,
             GetCharge(row*8+4,ibucket)/Munits::fC,
             GetCharge(row*8+5,ibucket)/Munits::fC,
             GetCharge(row*8+6,ibucket)/Munits::fC,
             GetCharge(row*8+7,ibucket)/Munits::fC,
             GetCharge(row*8+8,ibucket)/Munits::fC);
      
      printf("\n");
    }
  }
  printf("\n");
}


Float_t SimPmtM64::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.

  // Shortcut: saves a lot of unneccessary random number generation.
  // Note that there is no intrinsic charge resolution in the PMT
  // with this method.
  if(pe<=0) return 0;

  Float_t gain = fPixelGains[pixel];
  Float_t secemm = fPixelSecEmm[pixel];

  // 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);

  MSG("DetSim",Msg::kVerbose) << "Rolled: " << pe << " pe -> " << secondary_pe / secemm << " charge" << " with " << gain << " gain" << endl;

  // convert this into charge.
  Float_t q = (secondary_pe / secemm)  // charge in pe at anode
    * gain                        // Gain of pixel (unitless)
    * Munits::e_SI;               // e to charge
  
  return q;
}

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



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

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

    // Iterate over talking pixels.
    for(Int_t injPix = 1; injPix <= fNPixels; injPix++) {      
      SimPixelTimeBucket& pixtb_inj = pmttb.GetPixelBucket(injPix);
        
      // Only one spot (spot 1) per pixel.
      
      float injPE = pixtb_inj.GetPE(1);
      if( injPE > 0) { // There is charge in this spot.

          // Roll how many pe get crosstalked out of this pixel.

          float prob = fsCrosstalkTable->fOptTotProb[injPix];
          prob *= fsPmtScaleOpticalCrosstalk;         // Scale factor.
          float ttalk = fRandom->Poisson(prob*injPE); // Number of photons talked out
          int   ntalk = TMath::Nint(ttalk);           // ditto, as an integer

          for(int ipe=0;ipe<ntalk;ipe++) {
            // Move the photons around.
            
            // Choose the talked-to pixel.
            int whichPix = injPix;
            float r = fRandom->Uniform();
            for(Int_t xPix = 1; xPix <= fNPixels; xPix++) {      
              if(r<fsCrosstalkTable->fOptDist[injPix][xPix]) { 
                whichPix=xPix; 
                break;
                }
            }
            if(whichPix==injPix){
              MSG("DetSim",Msg::kWarning) << "Random number error choosing xtalk pixel.";
              continue;
            }
            
            // Debugging hack continued.
            //ofile << injPix << "\t" 
            //      << xpix << "\t" << pixtb_inj.GetPE(injSpot) << endl;
            
            SimPixelTimeBucket& pixtb = pmttb.GetPixelBucket(whichPix);

            // Move a photoelectron.
            MoveOpticalPE(1,  // move 1 pe
                          pixtb_inj,1,   // from pix, spot
                          pixtb,    1 ); // to pix, spot
            
            pixtb.SetTruthBit(DigiSignal::kCrosstalkOptical);

          } // Itr over talked pe 
         
        }
    } // Itr over inj pixel
     
    // Now we're done with all pixels in the bucket.
  } // it over buckets
  // whew!
  
  // Copy whatever didn't crosstalk.
  CopyPEtoPEXtalk();
}


Float_t SimPmtM64::GenElectricalCrosstalk( Int_t injPixel,
                                           Int_t xPixel, 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.

  // If no charge, do nothing:
  if(injCharge <= 0) return 0;

  float frac = fsCrosstalkTable->fElecFrac[injPixel][xPixel];
  float k    = fsCrosstalkTable->fElecK[injPixel][xPixel];

  frac *= fsPmtScaleChargeCrosstalk;

  // A fraction of the charge leaks.
  float charge = frac*injCharge;
  
  // The fraction might be variable, rather than constant.
  // If so, choose it from a distribution of width
  // proportional to the square root of the charge.

  // Don't bother with variability if it won't matter at the ~ 2 ADC level.
  const float kMinVariance2 = (2.0*Munits::fC)*(2.0*Munits::fC);

  float width2 = fabs(k*charge);
  if(width2 > kMinVariance2 ){
    float width = sqrt(width2); 
    return fRandom->Gaus(charge,width);
  };

  return charge;
}

---