Fiber.cc

Go to the documentation of this file.

#include <mdet/Fiber.h>

#include <mdet/Scintillator.h>
#include <mdet/MHierarchyInfo.h>
#include <mdet/Module.h>
//
#include <utl/Vector.h>
#include <utl/RandomEngine.h>
//
#include <fwk/RandomEngineRegistry.h>
//
#include <CLHEP/Random/RandPoisson.h>
#include <CLHEP/Random/RandExponential.h>
//
#include <cmath>

namespace {

  struct ScintShaper 
  {
    double GetWidth() const 
    {
      return fW;
    }
    double GetLength() const
    {
      return fL;
    }
    double GetHeight() const 
    {
      return fH;
    }
    void Box(const mdet::Scintillator&, double w, double l, double h) 
    {
      fW = w;
      fL = l;
      fH = h;
    }
    // No need to init them: the call to Box does that job.
    double fW;
    double fL;
    double fH;
  };

}

namespace mdet {

  const char* const 
  Fiber::kComponentName = MHierarchyInfo::kComponentsNames[4];
  
  const char* const  
  Fiber::kComponentId   = MHierarchyInfo::kComponentsIds[4];
 
  // Implement the pure virtual method from MPositionable.
  utl::CoordinateSystemPtr 
  Fiber::GetReferenceCoordinateSystem() 
    const
  {
    // We refer to the Scints's system.
    return GetModule().GetScintillatorFor(*this).GetLocalCoordinateSystem();
  }

  double 
  Fiber::GetOnManifoldLength() 
    const 
  {
    return GetData(fOnManifoldLength, "onManifoldLength"); 
  }

  double 
  Fiber::GetOnScintillatorLength() 
    const
  {
    // return GetData(fOnScintillatorLength, "onScintillatorLength"); 
    /*
     * Total overlap with its scintillator: it's glued in the scint.
     *
     * Note that this fact is not coupled to the actual position, which
     * is given by configuration: we leave those degrees of freedom, tough
     * assuming always this total overlap.
     */
    ScintShaper s;
    return GetModule().GetScintillatorFor(*this).VisitShape(s).GetLength();
  }

  double 
  Fiber::GetRadius() 
    const
  {
    return GetData(fRadius, "radius"); 
  }

  double 
  Fiber::GetNumericalAperture() 
    const
  {
    return GetData(fRadius, "numericalAperture"); 
  }

  double 
  Fiber::GetRefractionIndex() 
    const
  {
    return GetData(fRefractionIndex, "refractionIndex"); 
  }

  double 
  Fiber::GetDecayTime() 
    const
  {
    return GetData(fDecayTime, "decayTime"); 
  }

  double 
  Fiber::ComputeDecayDelay() 
    const
  {
    // This method gives sense to the bare number decay time: it's an exponential process.
        utl::RandomEngine& rnd = fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);

        double decayTime = CLHEP::RandExponential::shoot(& rnd.GetEngine(), GetDecayTime());

        return decayTime;
  }
  
  double 
  Fiber::GetDecayDelayMean() 
    const
  {
    // The mean is the parameter.
    return GetDecayTime();
  }

  double 
  Fiber::GetDecayDelayStdDev() 
    const
  {
    // For an exponential process the std-dev equals the mean value.
    return GetDecayDelayMean();
  }

  double 
  Fiber::GetAttenuationLengthA() 
    const
  {
    return GetData(fAttenuationLengthA, "attenuationLengthA"); 
  }
  
  double 
  Fiber::GetAttenuationLengthB() 
    const
  {
    return GetData(fAttenuationLengthB, "attenuationLengthB"); 
  }

  double 
  Fiber::GetAttenuationAmplitudeA() 
    const
  {
    return GetData(fAttenuationAmplitudeA, "attenuationAmplitudeA"); 
  }
  
  double 
  Fiber::GetAttenuationAmplitudeB() 
    const
  {
    return GetData(fAttenuationAmplitudeB, "attenuationAmplitudeB"); 
  }

  double
  Fiber::GetAttenuationLengthASiPM()
    const
  {
    return GetData(fAttenuationLengthASiPM, "attenuationLengthASiPM");
  }

  double
  Fiber::GetAttenuationLengthBSiPM()
    const
  {
    return GetData(fAttenuationLengthBSiPM, "attenuationLengthBSiPM");
  }

  double
  Fiber::GetAttenuationAmplitudeASiPM()
    const
  {
    return GetData(fAttenuationAmplitudeASiPM, "attenuationAmplitudeASiPM");
  }

  double
  Fiber::GetAttenuationAmplitudeBSiPM()
    const
  {
    return GetData(fAttenuationAmplitudeBSiPM, "attenuationAmplitudeBSiPM");
  }

  double 
  Fiber::GetAttenuationReferenceThickness() 
    const
  {
    return GetData(fAttenuationReferenceThickness, "attenuationReferenceThickness"); 
  }

  double
  Fiber::GetAttenuationReferenceEnergy()
    const 
  {
    return GetData(fAttenuationReferenceEnergy, "attenuationReferenceEnergy"); 
  }

  const char * const 
  Fiber::AttenuationReferenceTags[] = {
    "thickness", 
    "energy", 
  };

  Fiber::AttenuationReference 
  Fiber::GetAttenuationReference() 
    const
  {
    return AttenuationReferenceCreator::Create(GetData(fAttenuationReference, "attenuationReference"));
  }

  unsigned int 
  Fiber::ComputeSPENumber(double x, double l, double e) 
    const
  {

        double attenuationAmplitudeA;
        double attenuationAmplitudeB;
        double attenuationLengthA;
        double attenuationLengthB;

        if(!GetModule().IsSiPM()){
                attenuationAmplitudeA = GetAttenuationAmplitudeA();
                attenuationAmplitudeB = GetAttenuationAmplitudeB();
                attenuationLengthA = GetAttenuationLengthA();
                attenuationLengthB = GetAttenuationLengthB();

        }else{
                attenuationAmplitudeA = GetAttenuationAmplitudeASiPM();
                attenuationAmplitudeB = GetAttenuationAmplitudeBSiPM();
                attenuationLengthA = GetAttenuationLengthASiPM();
                attenuationLengthB = GetAttenuationLengthBSiPM();
        }

        double lambda0 = attenuationAmplitudeA * std::exp(-x / attenuationLengthA) +
                                                attenuationAmplitudeB * std::exp(-x / attenuationLengthB);

    double lambda = lambda0;
    switch (GetAttenuationReference()) {
      case eThickness:
        lambda *= l / GetAttenuationReferenceThickness();
      break;
      case eEnergy:
        lambda *= e / GetAttenuationReferenceEnergy();
      break;
    }
    utl::RandomEngine& rnd = fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);
    return CLHEP::RandPoisson::shoot(& rnd.GetEngine(), lambda);
  }


  Fiber::Fiber
    (int fId, const det::VManager::IndexMap& parentMap, const Module& parent) :
      MDetectorComponent<Fiber>::Type(fId, parentMap),
      det::MPositionable<Fiber>(*this),
      fModule(parent)
  {
    Register(fOnScintillatorLength);
    Register(fOnManifoldLength);
    Register(fNumericalAperture);
    Register(fRadius);
    Register(fRefractionIndex);
    Register(fDecayTime);
    Register(fAttenuationLengthA);
    Register(fAttenuationLengthB);
    Register(fAttenuationAmplitudeA);
    Register(fAttenuationAmplitudeB);
    Register(fAttenuationReferenceThickness);
  }

}