SimPmtM64Full.cxx

Go to the documentation of this file.

#include "SimPmtM64Full.h"
#include "SimPmtM64CrosstalkTable.h"
#include "SimPmtMaker.h"
#include "MessageService/MsgService.h"

CVSID("$Id: SimPmtM64Full.cxx,v 1.10 2007/03/01 17:22:11 rhatcher Exp $");

SimPmtMakerProxy<SimPmtM64Full> gSimPmtM64FullProxy("SimPmtM64Full");

ClassImp(SimPmtM64Full)

double SimPmtM64Full::fsChargeSlop = 1.6;

SimPmtM64Full::SimPmtM64Full(PlexPixelSpotId tube, 
            VldContext context,
            TRandom* random ) :
  SimPmtM64(tube,context,random)
{
  // Constructor.
}

Float_t SimPmtM64Full::PrevBucketFraction( Float_t dt )
{
  Float_t x = dt/Munits::ns;
  if(x > 3./fsChargeSlop) return 0;
  Float_t sx = x*fsChargeSlop;
  const Float_t k1 = -13.*exp(-3.);
  Float_t retval =  k1 + exp(-sx) * 
    ( 1 + sx * ( 1 + sx*(0.5 + sx/6.) ) );
  if((retval<0)||(retval>1)) 
    MSG("DetSim",Msg::kError) << "Prev Bucket leakage returned nonphysical value! Report this bug!" << endl;

  return retval;
}

Float_t SimPmtM64Full::NextBucketFraction( Float_t dt )
{
  // Calculates the fraction of charge 
  // from a single PE peak that will leak into the 
  // next PMT bucket, if the next pmt bucket is time (x) from
  // the time of the PE hit.

  // Assumes a functional form of t*t*t*exp(-t*1.6), with t in ns. 
  // See NUMI-L-757

  Float_t x = (dt/Munits::ns + 3./fsChargeSlop);

  Float_t sx = fsChargeSlop*x;
  Float_t retval = exp(-sx) * ( 6. + sx*(6. + sx*(3. + sx) ) ) / 6.;

  if((retval<0)||(retval>1)) 
    MSG("DetSim",Msg::kError) << "Prev Bucket leakage returned nonphysical value! Report this bug!" << endl;

  return retval;
}

void SimPmtM64Full::SimulateCharges()
{
  // This smushes charge into the next and previous QIE buckets.
  
  // It's assumed that digiPE occours at the peak of the pulse,
  // which has a shape of t^3 exp(-t*1.6).
  // Some of the charge gets put into the next and previous buckets.
  // Note that the truth info is lost: you will now get digits with 
  // no signal info... but that's the way it goes, I'm afraid.

  fTotalCharge = 0;
  // Iterate over all time buckets.
  for(SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    Int_t bucket_id = it.BucketId();
    SimPmtTimeBucket& bucket = it.Bucket();

    Float_t start_time = BucketToStartTime(bucket_id);
    Float_t stop_time  = BucketToStopTime(bucket_id);


    // Iterate over pixels.
    for(Int_t pix = 1; pix <= fNPixels; pix++) {      
      SimPixelTimeBucket& pixelBucket   = bucket.GetPixelBucket(pix);
      
      // Important: only look for next and prev buckets if there's 
      // something in here. The problem is that looking for the buckets
      // CREATES them, so if you always pull a ref to the next bucket,
      // the loop never exits.  But, if there is no work to do,
      // it doesn't create another bucket so we're safe.
      if(pixelBucket.GetTotalPEXtalk() >0) {

        SimPmtTimeBucket& nextBucket = GetBucket(bucket_id + 1);
        SimPixelTimeBucket& nextPixBucket = nextBucket.GetPixelBucket(pix);
        SimPmtTimeBucket& prevBucket = GetBucket(bucket_id - 1);
        SimPixelTimeBucket& prevPixBucket = prevBucket.GetPixelBucket(pix);

        for(Int_t spot = 1; spot <= fNSpots; spot++) {
          
          // Time to do it right.         
          // Loop over all the DigiPEs in this spot. Generate a charge for each.
          // Then put some of that charge in the previous or next time buckets,
          // depending on when the DigiPE hit relative to this bucket.
          
          SimPixelTimeBucket::PeList_t& pe_list = 
            pixelBucket.GetDigiPEXtalk(spot);   
          
          float npe  = pixelBucket.GetPEXtalk(spot);
          float qtot = GenChargeFromPE(pix,spot,npe);
          float q_per_pe = qtot/npe;

          SimPixelTimeBucket::PeList_t::iterator pe_itr;
          for(pe_itr = pe_list.begin(); pe_itr!=pe_list.end(); pe_itr++) {
            // Generate a charge.
            Float_t q = q_per_pe;
            
            // How much charge goes into the next bucket?
            Float_t frac_next =  NextBucketFraction( stop_time - (pe_itr->first) );
            Float_t frac_prev =  PrevBucketFraction( (pe_itr->first) - start_time );
            Float_t qnext = q* frac_next;
            Float_t qprev = q* frac_prev;

            // Add charge to next bucket, and set truth bits if there
            // is a non-trivial amount.
            nextPixBucket.AddCharge( qnext );
            if(qnext > 1.0*Munits::fC) 
              nextPixBucket.SetTruthBit(DigiSignal::kLeakFromPrevBucket);

            // Add charge to prev bucket, and set truth bits if there 
            // is a non-trivial amount.
            prevPixBucket.AddCharge( qprev );
            if(qprev > 1.0*Munits::fC) 
              prevPixBucket.SetTruthBit(DigiSignal::kLeakFromNextBucket);

            pixelBucket.AddCharge( q - qnext - qprev );
            bucket.AddTotalCharge(q);
            fTotalCharge += q;
          } // pes
        } // spots
      } // if
    } // pixels
  } // buckets.
}


void SimPmtM64Full::Config( Registry& config )
{
  double dtmp;
  if(config.Get("pmtChargeSlop",dtmp)) fsChargeSlop = dtmp; 
  SimPmtM64::Config(config); // call parent configuration.
}

---