MDetector/Pixel.cc

Go to the documentation of this file.

#include <mdet/Pixel.h>

#include <mdet/MHierarchyInfo.h>
//
#include <utl/RandomEngine.h>
#include <utl/MathConstants.h>
//
#include <fwk/RandomEngineRegistry.h>
//
#include <cmath>
//
#include <CLHEP/RandomObjects/RandMultiGauss.h>

namespace mdet {

  const char* const 
  Pixel::kComponentName = MHierarchyInfo::kComponentsNames[6];

  const char* const  
  Pixel::kComponentId   = MHierarchyInfo::kComponentsIds[6];

  const char * const 
  Pixel::PulseParametrizationTags[] = {
    "amplitudeAndCharge", 
    "amplitudeAndStandardDeviation", 
    "chargeAndStandardDeviation"
  };

  Pixel::Pixel
    (int pId, const det::VManager::IndexMap& parentMap, const PMT& parent) :
      MDetectorComponent<Pixel>::Type(pId, parentMap),
      fPMT(parent)
  { 
    Register(fCrossTalkNormalizationFactor);
    Register(fPulseParametrization);
    Register(fPulseAmplitudeMean);    
    Register(fPulseAmplitudeStdDev);
    Register(fPulseStdDevMean);
    Register(fPulseStdDevStdDev);
    Register(fPulseChargeMean);
    Register(fPulseChargeStdDev);
    Register(fPulseRelevantWidth);
  }

  const CrossTalkContainer&
  Pixel::GetNeighborsCrossTalkInside()
    const 
  {
    // Note that's shared: so once initialized no further loading
    // will be done. 
    return GetData<CrossTalkContainer,utl::ThrowOnZeroDereference,utl::ShadowPtr>(fNeighborsCrossTalk, "neighborsCrossTalkInside");
  }

  const CrossTalkContainer&
  Pixel::GetNeighborsCrossTalkSide()
    const 
  {
    // see the other getters.
    return GetData<CrossTalkContainer,utl::ThrowOnZeroDereference,utl::ShadowPtr>(fNeighborsCrossTalk, "neighborsCrossTalkSide");
  }

  const CrossTalkContainer&
  Pixel::GetNeighborsCrossTalkCorner()
    const 
  {
    // see the other getters.
    return GetData<CrossTalkContainer,utl::ThrowOnZeroDereference,utl::ShadowPtr>(fNeighborsCrossTalk, "neighborsCrossTalkCorner");
  }

  const IdsContainer&
  Pixel::GetChannelScintillatorFiberIds()
    const 
  {
    // see the other getters.
    return GetData<IdsContainer,utl::ThrowOnZeroDereference,utl::ShadowPtr>(fChannelScintillatorFiberIds, "channelScintillatorFiberIds");
  }

  Pixel::PulseParametrization 
  Pixel::GetPulseParametrization() 
    const
  {
    return PulseParametrizationCreator::Create(GetData(fPulseParametrization, "pulseParametrization"));
  }

  double 
  Pixel::GetPulseAmplitudeMean() 
    const
  {
    return GetData(fPulseAmplitudeMean, "pulseAmplitudeMean");
  }

  double 
  Pixel::GetPulseAmplitudeStdDev() 
    const 
  {
    return GetData(fPulseAmplitudeStdDev, "pulseAmplitudeStdDev");
  }
    
  double 
  Pixel::GetPulseStdDevMean() 
    const
  {
    return GetData(fPulseStdDevMean, "pulseStdDevMean");
  }

  double 
  Pixel::GetPulseStdDevStdDev() 
    const
  {
    return GetData(fPulseStdDevStdDev, "pulseStdDevStdDev");
  }

  double 
  Pixel::GetPulseChargeMean() 
    const
  {
    return GetData(fPulseChargeMean, "pulseChargeMean");
  }
  
  double 
  Pixel::GetPulseChargeStdDev() 
    const
  {
    return GetData(fPulseChargeStdDev, "pulseChargeStdDev");
  }

  double 
  Pixel::GetPulseParametersCorrelation() 
    const
  {
    return GetData(fPulseParametersCorrelation, "pulseParametersCorrelation");
  }

  double
  Pixel::GetPulseRelevantWidth()
    const
  {
    return GetData(fPulseRelevantWidth, "pulseRelevantWidth");
  }

  Pixel::SPE 
  Pixel::MakeSPEAt(double t)
    const
  {
    utl::RandomEngine& rnd = fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);
    // Provide spurious initialization for not-so-clever compilers: the following switch/case always enter in some case.
    double mu1    = 0;
    double mu2    = 0;
    double sigma1 = 0;
    double sigma2 = 0;
    PulseParametrization p = GetPulseParametrization();
    /*
     * Retrieve _the_ two parameters as of the chosen parametrization.
     * Take the name of the enumeration to denote the order of the parameters.
     */
    switch (p) {
      case eAmplitudeAndCharge:
        mu1    = GetPulseAmplitudeMean();
        sigma1 = GetPulseAmplitudeStdDev();
        mu2    = GetPulseChargeMean();
        sigma2 = GetPulseChargeStdDev();
      break;
      case eAmplitudeAndStandardDeviation:
        mu1    = GetPulseAmplitudeMean();
        sigma1 = GetPulseAmplitudeStdDev();
        mu2    = GetPulseStdDevMean();
        sigma2 = GetPulseStdDevStdDev();
      break;
      case eChargeAndStandardDeviation:
        mu1    = GetPulseChargeMean();
        sigma1 = GetPulseChargeStdDev();
        mu2    = GetPulseStdDevMean();
        sigma2 = GetPulseStdDevStdDev();
      break;
    }
    // Construct the parameters for the multivariate distro.
    CLHEP::HepSymMatrix covariance(2);
    covariance[0][0]  = sigma1 * sigma1;
    covariance[1][1]  = sigma2 * sigma2;
    // Recall that it's symmetric by construction: only one position set.
    covariance[1][0]  = GetPulseParametersCorrelation() * sigma1 * sigma2;
    // Vector with the two means.
    CLHEP::HepVector means(2);
    means[0] = mu1;
    means[1] = mu2;
    /*
     * Construct the object to be shoot several times.
     *  
     * Be careful to pass references to this class, if you pass a pointer
     * it accepts it anyway, but then the destructor delete the pointed-to
     * random engine: the common random engine gets deleted!
     */
    CLHEP::RandMultiGauss multiGauss(rnd.GetEngine(), means, covariance);
    /* 
     * Take random numbers.
     * But take care to have them positive: 
     * it's non sensical to have negative deviation
     * and, in this case, negative amplitude, or negative charge.
     */
    bool anyNegative;   // Loop go-on condition.
    CLHEP::HepVector r; // Collect the results.
    do {
      r = multiGauss.fire();
      anyNegative = false;
      for (int  i = 0; i < r.num_row() && !anyNegative; ++i) 
        anyNegative = r[i] < 0; // Note that we have a shortcut in the condition for the for-loop.
    } while (anyNegative);
    /*
     * Compute the parameters needed for the pulse, given the configured ones.
     * The parametrization is defined in terms of amplitude and standard-deviation;
     * if it's charge what's given, then obtain the final parameter according to the 
     * relation
     * Q = A * sigma * sqrt(2*pi)
     * where Q is the integral of the pulse all along the real axis.
     *
     * As before, provide spurious initialization to please compilers.
     */
    double a = 0;
    double s = 0;
    switch (p) {
      case eAmplitudeAndCharge:
        a = r[0];
        s = r[1] / r[0] / std::sqrt(utl::kTwoPi);
      break;
      case eAmplitudeAndStandardDeviation:
        a = r[0];
        s = r[1];
      break;
      case eChargeAndStandardDeviation:
        a = r[0] / r[1] / std::sqrt(utl::kTwoPi);
        s = r[1];
      break;
    }
    return SPE(*this, a, t, s);
  }

  void 
  Pixel::Update(bool invalidateData, bool invalidateComponents)
  {
    MDetectorComponent<Pixel>::Type::Update(invalidateData, invalidateComponents);
    // XXX The particular case of a shadow pointer doesn't fall into the VValidated
    // category...of course something could be abstracted...this is the only one case
    // so far....keep this here and that's it.
    if (invalidateData) {
      fNeighborsCrossTalk = 0;
      fChannelScintillatorFiberIds = 0;
    }
  }

}