#include <ArrayBM.h>

Inheritance diagram for ArrayBM:

Collaboration diagram for ArrayBM:

Public Member Functions
Constructor and Destructor
	ArrayBM (const std::string &name, ISvcLocator *svc)

virtual	~ArrayBM ()

AthService methods
virtual StatusCode	initialize () override final

IBeamIntensity methods
unsigned int	m_maxBunchCrossingPerOrbit
	max bunch crossings per orbit More...

Gaudi::Hive::ContextSpecificData< unsigned int >	m_t0Offset
	offset of the t0 wrto our intensity pattern More...

Gaudi::Property< std::uint64_t >	m_seed {this, "Seed", 0, "Seed for FastReseededPRNG. Zero seed switches to AthRNGSvc mode."}
	seed for FastReseededPRNG. Non-zero switches to using FastReseededPRNG. More...

Gaudi::Property< std::vector< float > >	m_intensityPatternProp
	user-defined intensity pattern More...

unsigned int	m_ipLength
	length of the intensity pattern More...

double *	m_intensityPattern
	normalized intensity pattern. C array to make clhep RandGeneral happy More...

CLHEP::RandGeneral *	m_biRandom
	shoot random number proportionally to m_intensityPattern More...

std::unique_ptr< std::discrete_distribution< unsigned int > >	m_t0Dist {nullptr}
	as with m_biRandom, but for FastReseededPRNG More...

ServiceHandle< IAthRNGSvc >	m_randomSvc {this, "RandomSvc", "AthRNGSvc","The random number service that will be used."}
	the service managing our random seeds/sequences More...

ATHRNG::RNGWrapper *m_rngWrapper	ATLAS_THREAD_SAFE {}

float	m_largestElementInPattern
	The largest value in the pattern assuming that the pattern has mean value 1.0. More...

int	m_emptyBunches
	Empty bunch option. More...

double *	m_signalPattern
	Additional array for keeping the locations we want signal in By default, will match the intensity pattern. More...

virtual float	normFactor (int iXing) const override final

virtual float	largestElementInPattern () const override final

virtual void	selectT0 (unsigned int run, unsigned long long event) override final

virtual unsigned int	getCurrentT0BunchCrossing () const override final

virtual unsigned int	getBeamPatternLength () const override final

Detailed Description

Definition at line 33 of file ArrayBM.h.

Constructor & Destructor Documentation

◆ ArrayBM()

ArrayBM::ArrayBM	(	const std::string &	name,
		ISvcLocator *	svc
	)

Definition at line 20 of file ArrayBM.cxx.

   : base_class(name,svc)
   , m_maxBunchCrossingPerOrbit(3564)
   , m_t0Offset(0)
   , m_intensityPatternProp()
   , m_ipLength(1)
   , m_intensityPattern(new double[m_ipLength])
   , m_biRandom(nullptr)
   , m_largestElementInPattern(1.0)
   , m_signalPattern(nullptr)
 {
   declareProperty("MaxBunchCrossingPerOrbit", m_maxBunchCrossingPerOrbit, "The number of slots in each LHC beam. Default: 3564.");
   declareProperty("IntensityPattern", m_intensityPatternProp,
                   "An array of floats containing the beam intensity distribution as a function of time in bins of 25ns. ArrayBM normalizes the distribution and uses it as a stencil to determine the relative intensity at each beam xing in the simulated range"
                   );
   declareProperty("EmptyBunchOption", m_emptyBunches=0, "Option for empty bunches.  0: normal treatment.  Positive N: first N BCIDs after filled.  Negative N: any empty BCID is allowed.");
   m_intensityPattern[0]=1.0;
 }

◆ ~ArrayBM()

ArrayBM::~ArrayBM ( )

virtual

Definition at line 39 of file ArrayBM.cxx.

 {
   delete [] m_intensityPattern;
   delete m_biRandom;
 }

Member Function Documentation

◆ getBeamPatternLength()

virtual unsigned int ArrayBM::getBeamPatternLength ( ) const

inlinefinaloverridevirtual

Definition at line 51 of file ArrayBM.h.

51 { return m_ipLength; }

◆ getCurrentT0BunchCrossing()

virtual unsigned int ArrayBM::getCurrentT0BunchCrossing ( ) const

inlinefinaloverridevirtual

Definition at line 50 of file ArrayBM.h.

50 { return m_t0Offset; }

◆ initialize()

StatusCode ArrayBM::initialize ( )

finaloverridevirtual

Definition at line 45 of file ArrayBM.cxx.

 {
   if (m_seed == 0u) {
     // Only need AthRNGSvc if seed is set to 0
     ATH_CHECK(m_randomSvc.retrieve());
     m_rngWrapper = m_randomSvc->getEngine(this);
   }
  
   // Need to copy to make modifications for empty bunches
   const std::vector<float>& rProp(m_intensityPatternProp.value());
   std::vector<float>::const_iterator pBegin(rProp.begin());
   std::vector<float>::const_iterator pEnd(rProp.end());
   m_ipLength = rProp.size();
   // Consistency checks
   if (m_ipLength>m_maxBunchCrossingPerOrbit)
     {
       ATH_MSG_ERROR("IntensityPattern length (" << m_ipLength << "), exceeds the maximum number of bunch crossings per orbit (" << m_maxBunchCrossingPerOrbit << ").");
       return StatusCode::FAILURE;
     }
  
   // Initialize the signal pattern if we need one different from the intensity pattern
   if (m_signalPattern) delete [] m_signalPattern;
   if (m_emptyBunches!=0){
     m_signalPattern = new double[m_ipLength];
   }
   // Modification for empty bunches option
   if (m_emptyBunches<0 || std::abs(m_emptyBunches)>static_cast<int>(m_ipLength)){
     // Easy case: Just flip all the bunches
     for (size_t i=0;i<m_ipLength;++i){
       if (rProp[i]>0.) m_signalPattern[i]=0.;
       else m_signalPattern[i]=1.;
     } // Loop over all bunches in the pattern
   } else if (m_emptyBunches>0){
     // Harder case: N BCIDs after filled
     int sinceFilled=0;
     // Set sinceFilled for the beginning of the pattern; assume we will not wrap (otherwise caught above)
     for (size_t i=m_ipLength-m_emptyBunches;i<m_ipLength;++i){
       if (rProp[i]>0) sinceFilled=0;
       else sinceFilled+=1;
     } // Done with the loop over previous BCIDs
     // Now do the loop setting the intensity pattern
     for (size_t i=0;i<m_ipLength;++i){
       if (rProp[i]>0){
         // Filled BCID.  Reset count, don't allow signal.
         sinceFilled=0;
         m_signalPattern[i]=0.;
       } else if (sinceFilled<m_emptyBunches){
         // First N BCIDs.  Increment count, allow signal.
         sinceFilled+=1;
         m_signalPattern[i]=1.;
       } else {
         // Beyond N BCIDs.  Increment count, don't allow signal.
         sinceFilled+=1;
         m_signalPattern[i]=0.;
       }
     } // Done with loop over previous BCIDs
   }
  
   // Normalise the pattern so that the non-zero elements average to 1.0
   float nonZeroElementCount(static_cast<float>(std::count_if(pBegin, pEnd, IsNonZero)));
   if(nonZeroElementCount<1.0)
     {
       ATH_MSG_ERROR("IntensityPattern has no non-zero elements!");
       return StatusCode::FAILURE;
     }
   float elementSum(static_cast<float>(std::accumulate(pBegin, pEnd,0.0)));
   float denominator(elementSum/nonZeroElementCount);
  
   // Normalise the pattern so that the highest element value is 1.0
   float maxElement(*(std::max_element(pBegin, pEnd)));
   float inv_maxElement = maxElement != 0 ? 1. / maxElement : 1;
  
   // Copy normalized intensity pattern from the property
   delete [] m_intensityPattern;
   m_intensityPattern = new double[m_ipLength];
   for (unsigned int i=0; i<m_ipLength; ++i)
     {
       if (rProp[i]<0.0)
         {
           ATH_MSG_ERROR("All IntensityPattern elements must be >=0. Please fix element #" << i );
           return StatusCode::FAILURE;
         }
       m_intensityPattern[i] = rProp[i] * inv_maxElement; // this ensures that the elements are all in the range [0,1]
     }
   // If we don't want to have signal in empty crossings, then our signal pattern is just the intensity pattern
   if (m_emptyBunches==0){
     m_signalPattern = m_intensityPattern;
   }
  
   // Will be used to convert values in the m_intensityPattern
   // from having max value 1.0 to having mean value 1.0
   if (denominator != 0) {
     m_largestElementInPattern = (maxElement/denominator);
   }
  
   // FIXME add a check that entry 0 is zero? In data, BCID=1 is always the
   // first filled bunch.
  
   if (m_seed == 0u) {
     delete m_biRandom;
     // the engine is created if not there already
     m_biRandom =
         new CLHEP::RandGeneral(m_signalPattern, m_ipLength,
                                /*IntType=*/1);  // discrete distribution
   } else {
     m_t0Dist =
         std::make_unique<std::discrete_distribution<unsigned int>>(
             m_intensityPattern, m_intensityPattern + m_ipLength);
   }
   return StatusCode::SUCCESS;
 }

◆ largestElementInPattern()

virtual float ArrayBM::largestElementInPattern ( ) const

inlinefinaloverridevirtual

Definition at line 48 of file ArrayBM.h.

48 { return m_largestElementInPattern; }

◆ normFactor()

float ArrayBM::normFactor ( int iXing ) const

finaloverridevirtual

Definition at line 179 of file ArrayBM.cxx.

 {
   unsigned int index = static_cast<unsigned int>((((iXing + static_cast<int>(m_t0Offset)) % static_cast<int>(m_ipLength)) + static_cast<int>(m_ipLength) )  % static_cast<int>(m_ipLength));
   //The array itself has max-value 1.0, but we want it to have mean value 1.0 (for filled bunch crossings), so we multiple by the m_largestElementInPattern.
   ATH_MSG_VERBOSE("normFactor for BCID " << iXing
                   << " (offset " << m_t0Offset
                   << " index " << index
                   << ") is = " <<  m_largestElementInPattern*m_intensityPattern[ index ]);
   return m_largestElementInPattern*m_intensityPattern[ index ];
 }

◆ selectT0()

void ArrayBM::selectT0	(	unsigned int	run,
		unsigned long long	event
	)

finaloverridevirtual

Definition at line 157 of file ArrayBM.cxx.

 {
   if (m_seed == 0u) {
     // Use AthRNGSvc
     const EventContext& ctx =
         Gaudi::Hive::currentContext();  // not ideal, but seems the cleanest
                                         // solution for now, as this call is
                                         // once per event.
     m_rngWrapper->setSeed("BEAMINT", ctx.slot(), event, run, ctx.evt());
     CLHEP::HepRandomEngine* rndmEngine = m_rngWrapper->getEngine(ctx);
     // m_biRandom->shoot() returns in range [0,1]
     m_t0Offset = static_cast<unsigned int>(
         floor((m_biRandom->shoot(rndmEngine) * m_ipLength) + 0.5));
   } else {
     FastReseededPRNG prng{m_seed.value(), run, event};
     m_t0Offset = (*m_t0Dist)(prng);
   }
   assert(m_intensityPattern[m_t0Offset % m_ipLength] > 0.0);  // just in case
   ATH_MSG_DEBUG("selectT0 offset for this event " << m_t0Offset);
   return;
 }

Member Data Documentation

◆ ATLAS_THREAD_SAFE

ATHRNG::RNGWrapper* m_rngWrapper ArrayBM::ATLAS_THREAD_SAFE {}

private

Definition at line 72 of file ArrayBM.h.

◆ m_biRandom

CLHEP::RandGeneral* ArrayBM::m_biRandom

private

shoot random number proportionally to m_intensityPattern

Definition at line 67 of file ArrayBM.h.

◆ m_emptyBunches

int ArrayBM::m_emptyBunches

private

Empty bunch option.

Default (0) means no special treatment of empty bunches, signal goes in filled bunches. Negative number means put signal in any empty bunch. Positive number means put signal in one of the first N bunches after a filled bunch.

Definition at line 82 of file ArrayBM.h.

◆ m_intensityPattern

double* ArrayBM::m_intensityPattern

private

normalized intensity pattern. C array to make clhep RandGeneral happy

Definition at line 65 of file ArrayBM.h.

◆ m_intensityPatternProp

Gaudi::Property<std::vector<float> > ArrayBM::m_intensityPatternProp

private

user-defined intensity pattern

Definition at line 61 of file ArrayBM.h.

◆ m_ipLength

unsigned int ArrayBM::m_ipLength

private

length of the intensity pattern

Definition at line 63 of file ArrayBM.h.

◆ m_largestElementInPattern

float ArrayBM::m_largestElementInPattern

private

The largest value in the pattern assuming that the pattern has mean value 1.0.

Multiplying by this converts values in the m_intensityPattern from having max value 1.0 to having mean value 1.0.

Definition at line 77 of file ArrayBM.h.

◆ m_maxBunchCrossingPerOrbit

unsigned int ArrayBM::m_maxBunchCrossingPerOrbit

private

max bunch crossings per orbit

Definition at line 55 of file ArrayBM.h.

◆ m_randomSvc

ServiceHandle<IAthRNGSvc> ArrayBM::m_randomSvc {this, "RandomSvc", "AthRNGSvc","The random number service that will be used."}

private

the service managing our random seeds/sequences

Definition at line 71 of file ArrayBM.h.

◆ m_seed

Gaudi::Property<std::uint64_t> ArrayBM::m_seed {this, "Seed", 0, "Seed for FastReseededPRNG. Zero seed switches to AthRNGSvc mode."}

private

seed for FastReseededPRNG. Non-zero switches to using FastReseededPRNG.

Definition at line 59 of file ArrayBM.h.

◆ m_signalPattern

double* ArrayBM::m_signalPattern

private

Additional array for keeping the locations we want signal in By default, will match the intensity pattern.

Definition at line 85 of file ArrayBM.h.

◆ m_t0Dist

std::unique_ptr<std::discrete_distribution<unsigned int> > ArrayBM::m_t0Dist {nullptr}

private

as with m_biRandom, but for FastReseededPRNG

Definition at line 69 of file ArrayBM.h.

◆ m_t0Offset

Gaudi::Hive::ContextSpecificData<unsigned int> ArrayBM::m_t0Offset

private

offset of the t0 wrto our intensity pattern

Definition at line 57 of file ArrayBM.h.

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

Public Member Functions

IBeamIntensity methods

Detailed Description

Constructor & Destructor Documentation

◆ ArrayBM()

◆ ~ArrayBM()

Member Function Documentation

◆ getBeamPatternLength()

◆ getCurrentT0BunchCrossing()

◆ initialize()

◆ largestElementInPattern()

◆ normFactor()

◆ selectT0()

Member Data Documentation

◆ ATLAS_THREAD_SAFE

◆ m_biRandom

◆ m_emptyBunches

◆ m_intensityPattern

◆ m_intensityPatternProp

◆ m_ipLength

◆ m_largestElementInPattern

◆ m_maxBunchCrossingPerOrbit

◆ m_randomSvc

◆ m_seed

◆ m_signalPattern

◆ m_t0Dist

◆ m_t0Offset