EventT0GeneratorModule.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
// Own header.
#include <generators/modules/EventT0GeneratorModule.h>
 
// framework aux
#include <framework/gearbox/Unit.h>
#include <framework/logging/Logger.h>
 
// root
#include <TRandom.h>
 
using namespace std;
using namespace Belle2;
 
//-----------------------------------------------------------------
//                 Register module
//-----------------------------------------------------------------
 
REG_MODULE(EventT0Generator);
 
//-----------------------------------------------------------------
//                 Implementation
//-----------------------------------------------------------------
 
EventT0GeneratorModule::EventT0GeneratorModule() : Module()
 
{
  // set module description
  setDescription("Module generates discrete event t0 in ~4ns steps (bunch spacing) "
                 "according to double gaussian distribution and adds it to the "
                 "production and decay times of MCParticles. This means that after "
                 "this module the time origin (t = 0) is set to what L1 trigger "
                 "would give as the collision time. In case of cosmics, the L1 trigger"
                 "jitter is generated according to a continuos double gaussian distribution");
 
  setPropertyFlags(c_ParallelProcessingCertified);
 
  // Add parameters
  addParam("coreGaussWidth", m_coreGaussWidth, "sigma of core gaussian [ns]", double(5.6));
  addParam("tailGaussWidth", m_tailGaussWidth, "sigma of tail gaussian [ns]", double(14.5));
  addParam("tailGaussFraction", m_tailGaussFraction,
           "fraction (by area) of tail gaussian", double(0.08));
  addParam("fixedT0",  m_fixedT0,
           "If set, a fixed event t0 is used instead of simulating the bunch timing.", m_fixedT0);
  addParam("maximumT0",  m_maximumT0,
           "If set, randomize between -maximum and maximum.",
           m_maximumT0);
  addParam("isCosmics", m_isCosmics,
           "if True simulate L1 jitter for cosmics", bool(false));
  addParam("coreGaussWidthCosmics", m_coreGaussWidthCosmics, "sigma of core gaussian for cosmics [ns]", double(13));
  addParam("tailGaussMeanCosmics", m_tailGaussMeanCosmics, "mean of tail gaussian for cosmics [ns]", double(28));
  addParam("tailGaussWidthCosmics", m_tailGaussWidthCosmics, "sigma of tail gaussian for cosmics [ns]", double(15));
  addParam("tailGaussFractionCosmics", m_tailGaussFractionCosmics,
           "fraction (by area) of tail gaussian for cosmics", double(0.09));
 
}
 
 
void EventT0GeneratorModule::initialize()
{
  m_mcParticles.isRequired();
  m_initialParticles.registerInDataStore();
  m_simClockState.registerInDataStore();
 
  if (not std::isnan(m_maximumT0) and not std::isnan(m_fixedT0)) {
    B2ERROR("You can not set both the maximum T0 and the fixed T0 option.");
  }
}
 
 
void EventT0GeneratorModule::event()
{
 
 
  double eventTime = 0;
  int bucket = 0;
  int beamRevo9Cycle = 0;
  int relBucketNo = 0;
  //if m_isCosmics we need to randomize revo9count,
  //at least for the SVD trigger bin
  //if !m_isCosmics the revo9Count value is overwritten
  unsigned revo9range = (m_bunchStructure->getRFBucketsPerRevolution() / 4) * 9;
  int revo9count = gRandom->Integer(revo9range);
 
  if (!m_isCosmics) {
    // generate bucket number w.r.t revo9 marker
 
    bucket = m_bunchStructure->generateBucketNumber(); // [RF clock]
    beamRevo9Cycle = gRandom->Integer(9); // number of beam revolutions w.r.t revo9 marker
    int bucketRevo9 = bucket + beamRevo9Cycle * m_bunchStructure->getRFBucketsPerRevolution(); // [RF clock]
 
    // generate L1 time jitter
 
    double timeJitter = 0;
    if (not std::isnan(m_maximumT0)) {
      timeJitter = -m_maximumT0 + (2 * m_maximumT0) * gRandom->Rndm();
    } else if (not std::isnan(m_fixedT0)) {
      timeJitter = m_fixedT0;
    } else {
      double sigma = m_coreGaussWidth;
      if (gRandom->Rndm() < m_tailGaussFraction) sigma = m_tailGaussWidth;
      timeJitter = gRandom->Gaus(0., sigma);
    }
 
    // calculate revo9count (system clock ticks since revo9 marker as determined by L1 trigger)
 
    double bucketTimeSep = 1 / m_clockSettings->getAcceleratorRF();
    relBucketNo = round(timeJitter / bucketTimeSep); // RF clock
    revo9count = (bucketRevo9 + relBucketNo) / 4; // system clock
 
    // calculate collision time w.r.t L1 trigger
 
    eventTime = (bucketRevo9 - revo9count * 4) * bucketTimeSep;
 
  } else {
 
    //compute jitter for cosmics, .i.e. no filling pattern
    double sigma = m_coreGaussWidthCosmics;
    if (gRandom->Rndm() < m_tailGaussFractionCosmics) sigma = m_tailGaussWidthCosmics;
    eventTime = gRandom->Gaus(0., sigma);
  }
 
  // correct MC particle times according to generated collision time
 
  for (auto& particle : m_mcParticles) {
    particle.setProductionTime(particle.getProductionTime() + eventTime);
    particle.setDecayTime(particle.getDecayTime() + eventTime);
  }
 
  // store collision time to MC initial particles (e.g. beam particles)
 
  if (not m_initialParticles.isValid()) m_initialParticles.create();
  m_initialParticles->setTime(eventTime);
 
  // store revo9count (modulo range) in order to be distibuted to sub-detectors
  revo9count %= revo9range;
  if (revo9count < 0) revo9count += revo9range;
 
  m_simClockState.create();
  m_simClockState->setRevo9Count(revo9count);
  m_simClockState->setBucketNumber(bucket);
  m_simClockState->setBeamCycleNumber(beamRevo9Cycle);
  m_simClockState->setRelativeBucketNo(relBucketNo);
}