TileDigitsFromPulse.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
//*****************************************************************************
//  Filename : TileDigitsFromPulse.cxx
//  Author   : Simon Molander
//  Created  : February 2013
//
//  DESCRIPTION
// 
//  Create TileDigits from simulated pulses. 
//
//  HISTORY:
//
//  BUGS:
//
//*****************************************************************************
 
// Tile includes
#include "TileSimAlgs/TileDigitsFromPulse.h"
#include "TileEvent/TileDigits.h"
#include "TileEvent/TileMutableDigitsContainer.h"
#include "TileEvent/TileMutableRawChannelContainer.h"
#include "TileIdentifier/TileHWID.h"
#include "TileConditions/TileInfo.h"
#include "TileCalibBlobObjs/TileCalibUtils.h"
 
//Simulator includes
#include "TilePulseSimulator/TileSampleGenerator.h"
#include "TilePulseSimulator/TileSampleBuffer.h"
#include "TilePulseSimulator/TilePulseShape.h"
 
// Athena includes
#include "AthAllocators/DataPool.h"
#include "PathResolver/PathResolver.h"
//Random number service
#include "AthenaKernel/IAthRNGSvc.h"
#include "AthenaKernel/RNGWrapper.h"
 
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Units/SystemOfUnits.h>
 
//Root includes
#include "TRandom3.h"
#include "TFile.h"
#include "TH1F.h"
#include "TKey.h"
#include "TF1.h"
 
#include <cstdlib>
 
//C++ STL includes
#include <vector>
using CLHEP::RandGaussQ;
using CLHEP::RandFlat;
 
//
// Constructor
//
TileDigitsFromPulse::TileDigitsFromPulse(const std::string& name, ISvcLocator* pSvcLocator) :
  AthAlgorithm(name, pSvcLocator)
{
 
    declareProperty("ImperfectionMean", m_imperfectionMean = 1.01);
    declareProperty("ImperfectionRms", m_imperfectionRms = 0.02);
    declareProperty("InTimeAmp", m_inTimeAmp = 300.);
    declareProperty("OutOfTimeAmp", m_ootAmp = 150.);
    declareProperty("InTimeOffset", m_itOffset = 0.);
    declareProperty("OutOfTimeOffset", m_ootOffset = 50.);
    declareProperty("OutOfTimeOffsetHistogramFile", m_ootOffsetFileName = "");
    declareProperty("OutOfTimeOffsetHistogramName", m_ootOffsetHistName = "");
    declareProperty("UseGaussNoise", m_gaussNoise = kFALSE);
    declareProperty("GaussNoiseAmpOne", m_GNAmpOne = 1 / 1.039);
    declareProperty("GaussNoiseSigmaOne", m_GNSigmaOne = 1.6);
    declareProperty("GaussNoiseAmpTwo", m_GNAmpTwo = 0.039);
    declareProperty("GaussNoiseSigmaTwo", m_GNSigmaTwo = 3.6);
    declareProperty("UseInTimeAmpDist", m_useItADist = kFALSE);
    declareProperty("UseOutOfTimeAmpDist", m_useOotADist = kFALSE);
    declareProperty("InTimeAmpDistFileName", m_itADistFileName = "");
    declareProperty("InTimeAmpPulseProb", m_itAPulseProb = 0);
    declareProperty("OutOfTimeAmpDistFileName", m_ootADistFileName = "");
    declareProperty("PileUpFraction", m_pileUpFraction = 1);
    declareProperty("GaussianC2CPhaseVariation", m_gausC2C = 0);
    declareProperty("ChannelSpecificPedestal", m_chanPed = kFALSE);
    declareProperty("ChannelSpecificNoise", m_chanNoise = kFALSE);
    declareProperty("PedestalValueHG", m_ped_HG = 100);
    declareProperty("PedestalValueLG", m_ped_LG = 100);
    declareProperty("AmpDistLowerLimit", m_AmpDistLowLim = 135);
    declareProperty("InTimeAmpDistHistogramName", m_itADistHistName = "h_Eopt_hi");
    declareProperty("OutOfTimeAmpDistHistogramName", m_ootADistHistName = "h_Eopt_hi");
 
    declareProperty("SimulatePileUpWithPoiss", m_simPUwPoisson = kFALSE); 
    declareProperty("AvgMuForPileUpSimulation", m_avgMuForPU = 40);
    declareProperty("PileUpAmpDistFileName", m_pileupAmpDistFileName = "");
 
    declareProperty("RandomSeed", m_seed = 4357);
    declareProperty("BunchSpacing", m_BunchSpacing = 25.); // 25, 50 or 75
    declareProperty("SimulateQIE", m_simQIE = kFALSE);
    declareProperty("SimulatePulseChain",m_simPulseChain = kFALSE);
 
    declareProperty("TileInfoName", m_infoName = "TileInfo");
    declareProperty("TilePhaseII", m_PhaseII = kFALSE);
    declareProperty("Bigain", m_bigain = kFALSE);
    declareProperty("NSamples", m_nSamples = 7); 
    declareProperty("NPulses", m_nPul = 21);
 
    //Initialisations
    m_ps[0] = new TilePulseShape(msgSvc(), "TilePulseShapeLo"); //Low Gain
    m_ps[1] = new TilePulseShape(msgSvc(), "TilePulseShapeHi"); //High Gain
 
    
    m_itFile = new TFile();
    m_itDist = new TH1F();
    m_ootFile = new TFile();
    m_ootDist = new TH1F();
    m_ootOffsetDist = new TH1F();
    m_ootOffsetFile = new TFile();
 
    m_pileup_AmpDistFile = new TFile();
 
    m_useOffsetHisto = kFALSE;
 
}
 
TileDigitsFromPulse::~TileDigitsFromPulse() {
 
    delete m_ootOffsetFile;
    delete m_ootFile;
    delete m_itFile;
    delete m_pileup_AmpDistFile;
    
    delete m_ps[0];
    delete m_ps[1];
}
 
//
// Alg standard initialize function
//
StatusCode TileDigitsFromPulse::initialize() {
 
    ATH_MSG_DEBUG("in initialize()");
 
    m_buf = new TileSampleBuffer(m_nSamples, -(m_nSamples-1)*25/2, 25.);
    
    m_tsg = new TileSampleGenerator(m_ps[0], m_buf, false); //Set third parameter to true for debug of the sum of pulses
 
    m_nPul_eff = (m_nPul - 1) / 2; //Used for symetrization of PU in computation
 
    ATH_CHECK(detStore()->retrieve(m_tileHWID, "TileHWID"));
    ATH_CHECK(detStore()->retrieve(m_tileInfo, m_infoName));
    m_i_ADCmax = m_tileInfo->ADCmax();
    ATH_MSG_DEBUG("Max ADC counts TileDigitsFromPulse: " << m_i_ADCmax);
 
        ATH_CHECK(m_tileToolNoiseSample.retrieve());
 
        ATH_CHECK(m_digitsContainerKey.initialize());
        ATH_MSG_INFO("Output digits container: " <<  m_digitsContainerKey.key());
 
        ATH_CHECK(m_rawChannelContainerKey.initialize());
    ATH_MSG_INFO("Output raw channel container: " << m_rawChannelContainerKey.key());
 
 
    //Build pulse shapes
    m_ps[0]->setPulseShape(m_tileInfo->digitsFullShapeLo());
    m_ps[1]->setPulseShape(m_tileInfo->digitsFullShapeHi());
 
    m_random = std::make_unique<TRandom3>(m_seed);
 
    //Initialise distribution histograms if in use
    if (m_useItADist) {
        if (m_itADistFileName.size() == 0) {
            m_itADistFileName = PathResolver::find_file("Distributions_small_h2000_177531_JetTauEtmiss.root", "DATAPATH");
            if (m_itADistFileName.size() == 0) {
                ATH_MSG_FATAL("Could not find input file Distributions_small_h2000_177531_JetTauEtmiss.root");
                return StatusCode::FAILURE;
            }
        }
        if (makeDist(m_itFile, m_itDist, m_itADistFileName, m_itADistHistName) == kFALSE)
            return StatusCode::FAILURE;
        ATH_MSG_DEBUG("Made in-time distribution");
    } else
        delete m_itDist;
 
    if (m_simPUwPoisson){
 
        if (m_pileupAmpDistFileName.size() == 0) {
            m_pileupAmpDistFileName = PathResolver::find_file("Distributions_MB_minbias_inelastic_lowjetphoton_e8314_e7400_s3508.root", "DATAPATH");
            if (m_pileupAmpDistFileName.size() == 0 ) {
                ATH_MSG_FATAL("Could not find input file Distributions_MB_minbias_inelastic_lowjetphoton_e8314_e7400_s3508.root");
                return StatusCode::FAILURE;
            }
 
        }
        if (makeDist(m_pileup_AmpDistFile, m_pileup_AmpDists, m_pileupAmpDistFileName) == kFALSE)
            return StatusCode::FAILURE;
        ATH_MSG_DEBUG("Made PU amp distributions for each partition and channel");
 
    }
 
    if (m_useOotADist) {
        if (m_ootADistFileName.size() == 0) {
            m_ootADistFileName = PathResolver::find_file("Distributions_small_h2000_177531_ZeroBias.root", "DATAPATH");
            if (m_ootADistFileName.size() == 0) {
                ATH_MSG_FATAL("Could not find input file Distributions_small_h2000_177531_ZeroBias.root");
                return StatusCode::FAILURE;
            }
        }
        if (makeDist(m_ootFile, m_ootDist, m_ootADistFileName, m_ootADistHistName) == kFALSE)
            return StatusCode::FAILURE;
        ATH_MSG_DEBUG("Made Oot distribution");
    } else
        delete m_ootDist;
 
    //Initialise timing offset distribution. If filename is empty, use static offset
    if (m_ootOffsetFileName.size() != 0) {
        m_ootOffsetFile = TFile::Open(m_ootOffsetFileName.c_str());
        if (m_ootOffsetFile->IsZombie()) {
            ATH_MSG_WARNING("Error reading offset timing distribution from " << m_ootOffsetFileName << ". Using static timing offset.");
        } else {
            TKey *key = m_ootOffsetFile->FindKey(m_ootOffsetHistName.c_str());
            if (key == 0) {
                ATH_MSG_WARNING("Histogram " << m_ootOffsetHistName << " not found in file " << m_ootOffsetFileName << ". Using static timing offset.");
            } else {
                m_ootOffsetDist = (TH1F*) m_ootOffsetFile->Get(m_ootOffsetHistName.c_str());
                m_useOffsetHisto = kTRUE;
            }
        }
    }
 
    //Start the random number service used to create channel specific noise
    if (!m_rndmSvc.retrieve().isSuccess()) {
        ATH_MSG_FATAL("Could not initialize find Random Number Service.");
        return StatusCode::FAILURE;
    }
    if (m_chanNoise)
        m_gaussNoise = kFALSE; //Make sure channel noise overrides gaussian noise.
 
    m_PUAmp.clear();
    m_PUAmp.resize(4);
    for (int ros = 1; ros < 5; ++ros){
      m_PUAmp[ros-1].clear();
      m_PUAmp[ros-1].resize(64);
 
      for (int drawer = 0; drawer < 64; ++drawer){
        m_PUAmp[ros-1][drawer].clear();
        m_PUAmp[ros-1][drawer].resize(48);
        for (int channel = 0; channel < 48; channel++){
          m_PUAmp[ros-1][drawer][channel].clear();
          m_PUAmp[ros-1][drawer][channel].resize(m_nPul);
        }
      }
    }
 
    for (int ros = 1; ros < 5; ++ros) { // initialize the vector of PU amplitudes
      for (int drawer = 0; drawer < 64; ++drawer){
        for (int channel = 0; channel < 48; ++channel) {
          addPileUp(m_inTimeAmp, 1, ros, drawer, channel); // Initialized for HG, for LG  use the same divided by the HG/LG ratio
        }
      }
    }
 
    ATH_MSG_DEBUG("initialize() successful");
 
    return StatusCode::SUCCESS;
}
/*==========================================================================*/
//
// Begin Execution Phase.
//
StatusCode TileDigitsFromPulse::execute() {
 
    ATH_MSG_DEBUG("in execute()");
 
    const EventContext& ctx = Gaudi::Hive::currentContext();
 
    // Prepare RNG service
    ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this, m_randomStreamName);
    rngWrapper->setSeed( m_randomStreamName, ctx );
 
    // Create new container for digits
    auto digitsContainer = std::make_unique<TileMutableDigitsContainer>(true,
                                                                            TileFragHash::Digitizer,
                                                                            TileRawChannelUnit::ADCcounts,
                                                                            SG::VIEW_ELEMENTS);
 
        ATH_CHECK( digitsContainer->status() );
 
    //Create RawChannel for truth values.
    auto rawChannelContainer = std::make_unique<TileMutableRawChannelContainer>(true, m_rChType, m_rChUnit);
        ATH_CHECK( rawChannelContainer->status() );
 
    DataPool < TileDigits > tileDigitsPool(m_tileHWID->adc_hash_max());
 
    double tFit = 0, ped = 100; //Settings for simulation 
 
    TF1 *pdf = new TF1();
    TF1 *pdf_PhaseI = new TF1();
    TF1 *pdf_lo = new TF1();
    TF1 *pdf_hi = new TF1();
    if (!m_simQIE) {
        
 
        if(m_PhaseII){
          Double_t sigma_lo = 1; //Noise value obtained from ATL-COM-TILECAL-2020-031
            pdf_lo = new TF1("pdf_lo","(1/(sqrt(2*pi)*[0])) * (exp(-0.5*(x/[0])**2)/(sqrt(2*pi)*[0]))", -100, 100);
            pdf_lo->SetParameter(0,sigma_lo);
 
            Double_t sigma_hi = 2.5; //Noise value obtained from ATL-COM-TILECAL-2020-031
            pdf_hi = new TF1("pdf_hi","(1/(sqrt(2*pi)*[0])) * (exp(-0.5*(x/[0])**2)/(sqrt(2*pi)*[0]))", -100, 100);
            pdf_hi->SetParameter(0,sigma_hi);
        }
        else{
            //Noise pdf for general noise. Maybe use as a member and put in init.
            //pdf = new TF1("pdf", "[0] * (Gaus(x,0,[1]) + [2] * Gaus(x,0,[3]))", -100., 100.); //Root goes not like "Gaus"     
 
            pdf_PhaseI = new TF1("pdf_PhaseI", "[0] * (exp(-0.5*(x/[1])**2)/(sqrt(2*pi)*[1]) + [2] *exp(-0.5*(x/[3])**2)/(sqrt(2*pi)*[3]))", -100., 100.);
            pdf_PhaseI->SetParameters(m_GNAmpOne, m_GNSigmaOne, m_GNAmpTwo, m_GNSigmaTwo);
        }
    }
 
    std::vector<float> samples(m_nSamples);
 
    double Rndm[16]; // Can't use variable size array,
    double Rndm_dG[1]; // uniform random number for the double gaussian
 
    ATH_MSG_DEBUG("Starting loop");
    int gain = 1;
    double n_inTimeAmp = 0.0; 
    float sample = 0.0;
 
    for (int ros = 1; ros < 5; ++ros) {
        for (int drawer = 0; drawer < 64; ++drawer) {
            unsigned int drawerIdx = TileCalibUtils::getDrawerIdx(ros, drawer);
            for (int channel = 0; channel < 48; ++channel) {
 
                if (!m_simQIE && !m_simPulseChain) { //3-in-1 or FENICS is simulated below
 
                    if(m_PhaseII){
                        ATH_MSG_VERBOSE("executing FENICS code");
                    }
                    else{
                        ATH_MSG_VERBOSE("executing 3-in-1 code");
                    }
                    
 
                    bool isHGSaturated = false;
 
                    for (int igain = 1; igain >= 0; igain--) {
 
                        gain = igain;
                        if (m_chanPed){
                          ped = m_tileToolNoiseSample->getPed(drawerIdx, channel, gain, TileRawChannelUnit::ADCcounts, ctx);
                        }
                        else{
                          ped = (gain == 1) ? m_ped_HG : m_ped_LG;
                        }
 
                        if (gain == 1) {
                            n_inTimeAmp = m_useItADist ? m_itDist->GetRandom() : m_inTimeAmp;
 
                            if (m_random->Rndm() >= m_pileUpFraction){
                                    m_ootAmp = 0; //Set oot amplitude to 0 if no pile-up.
                            }
                            tFit = m_random->Gaus(0., m_gausC2C); //C2C phase variation
                            double deformatedTime = m_random->Gaus(m_imperfectionMean, m_imperfectionRms); //Widening of pulseshape
                            m_ps[gain]->scalePulse(deformatedTime, deformatedTime); // Deformation of pulse shape by changing its width
                            //if(m_useOffsetHisto) m_ootOffset = m_ootOffsetDist->GetRandom();  //OLD Remove for 7 samples -> BunchSpacing
 
                            // Make sure m_PUAmp[ros-1][drawer][channel] has m_nPul elements, all zero
                            m_PUAmp[ros-1][drawer][channel].assign(m_nPul, 0.0);
                            
                            // Fill amplitudes (including pile-up) for HG
                            addPileUp(n_inTimeAmp, gain, ros, drawer, channel);
 
                        } else {
                            double deformatedTime = m_random->Gaus(m_imperfectionMean, m_imperfectionRms); //Widening of pulseshape
                            m_ps[gain]->scalePulse(deformatedTime, deformatedTime); // Deformation of pulse shape by changing its width
 
                            double scaleFactor = (m_PhaseII ? 40.0 : 64.0);
                            n_inTimeAmp /= scaleFactor;
                            for (auto &ampValue : m_PUAmp[ros-1][drawer][channel]) {
                                    ampValue /= scaleFactor;
                            }
                        }
 
                        if(m_PhaseII && m_gaussNoise){
                            pdf = (gain==1) ? pdf_hi : pdf_lo;
                        }
                        else if(m_gaussNoise){
                            pdf = pdf_PhaseI;
                        }
 
                        m_tsg->setPulseShape(m_ps[gain]);
                        m_tsg->fillNSamples(tFit, ped, n_inTimeAmp, m_PUAmp[ros-1][drawer][channel], pdf, m_gaussNoise, m_itOffset, m_nSamples, m_nPul); // Sum of Intime + PU pulses
 
                        samples.clear();
                        samples.resize(m_nSamples);
                        m_buf->getValueVector(samples);
 
                        if (m_chanNoise) {
                            double Hfn1 = m_tileToolNoiseSample->getHfn1(drawerIdx, channel, gain, ctx);
                            double Hfn2 = m_tileToolNoiseSample->getHfn2(drawerIdx, channel, gain, ctx);
                            double Norm = m_tileToolNoiseSample->getHfnNorm(drawerIdx, channel, gain, ctx);
                            RandGaussQ::shootArray(*rngWrapper, samples.size(), Rndm, 0.0, 1.0);
                            RandFlat::shootArray(*rngWrapper, 1, Rndm_dG, 0.0, 1.0);
                            for (unsigned int js = 0; js < samples.size(); ++js) {
                                //using the same gaussian(sigma) for all samples in one channel in one event
                                if (Rndm_dG[0] < Norm)
                                    samples[js] += (float) Hfn1 * Rndm[js];
                                else
                                    samples[js] += (float) Hfn2 * Rndm[js];
                            }
                        }
 
                        for (unsigned int i = 0; i < samples.size(); ++i) {
                          if (samples[i] >= m_i_ADCmax){
                            isHGSaturated = true;
                            if(m_bigain) samples[i]=m_i_ADCmax;
                          }
                            
                        }
                        
                        if(!m_bigain && !isHGSaturated){
                                break;
                        }
                        
                        ATH_MSG_VERBOSE("New ADC " << ros << "/" << drawer << "/" << channel << "/   saving gain  " << gain);
 
                        TileDigits * digit = tileDigitsPool.nextElementPtr();
                        *digit = TileDigits (m_tileHWID->adc_id(ros, drawer, channel, gain),
                                     std::move(samples));
 
                        ATH_CHECK( digitsContainer->push_back(digit) );
                          
                        auto rawChannel = std::make_unique<TileRawChannel>(digit->adc_HWID(),
                                                     n_inTimeAmp,
                                                     tFit,
                                                     m_ootAmp,
                                                     m_ootOffset);
                          
                        ATH_CHECK( rawChannelContainer->push_back(std::move(rawChannel)) );
 
                    }
 
                    if(!m_bigain){
                            ATH_MSG_VERBOSE("New ADC " << ros << "/" << drawer << "/" << channel << "/   saving gain  " << gain);
 
                        TileDigits * digit = tileDigitsPool.nextElementPtr();
                        *digit = TileDigits (m_tileHWID->adc_id(ros, drawer, channel, gain),
                                     std::move(samples));
 
                        ATH_CHECK( digitsContainer->push_back(digit) );
 
                        auto rawChannel = std::make_unique<TileRawChannel>(digit->adc_HWID(),
                                                     n_inTimeAmp,
                                                     tFit,
                                                     m_ootAmp,
                                                     m_ootOffset);
 
                        ATH_CHECK( rawChannelContainer->push_back(std::move(rawChannel)) );
                    }
 
                } else if (m_simQIE) { //QIE is simulated here --------------------------------------------
 
                    //ATH_MSG_DEBUG("executing QIE code");
 
                    gain = 1; //This is just a place holder. The gain is not used in QIE.
                    n_inTimeAmp = m_useItADist ? m_itDist->GetRandom() : m_inTimeAmp;
                    //if (random->Rndm() >= m_pileUpFraction) //m_pileUpFraction is 1 by default
                    m_ootAmp = 0; //Set oot amplitude to 0 if no pile-up.
                    tFit = 0; //TODO: Introduce jitter of the PMT pulse; random->Gaus(0., m_gausC2C); //C2C phase variation
 
                    //Pileup samples
                    //m_PUAmp.clear();
                    //m_PUAmp.resize(nPul);
                    float my_PUAmp[7] = {0}; //I use an array to store the energies/charges of the out-of-time pulses
 
                    for (int i = 0; i < 7; i++)
                        if ((((i - 3) * 25) % (int) m_BunchSpacing) == 0) {
                            if (i != 3) { //index 3 corresponds to the in-time pulse, the signal
                                my_PUAmp[i] = m_useOotADist ? m_ootDist->GetRandom() : m_ootAmp; //out-of-time pulses
                            } else {
                                my_PUAmp[i] = 0;
                            }
                        }
 
                    //fill7SamplesQIE(float t0, float amp_it, float *amp_pu, bool addNoise);
                    m_tsg->fill7SamplesQIE((float) n_inTimeAmp, my_PUAmp); // Sum of In time + out-of-time PU pulses
 
                    samples.clear();
                    samples.resize(m_nSamples);
                    m_buf->getValueVector(samples);
 
                    ATH_MSG_VERBOSE("New ADC " << ros << "/" << drawer << "/" << channel << "/   saving gain  " << gain);
 
                    TileDigits * digit = tileDigitsPool.nextElementPtr();
                    *digit = TileDigits (m_tileHWID->adc_id(ros, drawer, channel, gain),
                             std::move(samples));
                
                    ATH_CHECK( digitsContainer->push_back(digit) );
                
                    auto rawChannel = std::make_unique<TileRawChannel>(digit->adc_HWID(),
                                           n_inTimeAmp,
                                           tFit,
                                           m_ootAmp,
                                           m_ootOffset);
                
                    ATH_CHECK( rawChannelContainer->push_back(std::move(rawChannel)) );
                }
                else if (m_simPulseChain) {
                        ATH_MSG_VERBOSE("executing chain-of-pulses code");
 
                    bool isHGSaturated = false;
 
                    for (int igain = 1; igain >= 0; --igain) {
                              gain = igain;
                              n_inTimeAmp = 0.0;
                              m_ootAmp = 0.0;
                              m_ootOffset = 0.0;
 
                          if (m_random->Rndm() < m_itAPulseProb){
                            n_inTimeAmp = m_useItADist ? m_itDist->GetRandom() : m_inTimeAmp;
                          } else{
                            n_inTimeAmp = 0;
                          }
 
                          // PDF logic for noise
                          if (m_gaussNoise) {
                            if (m_PhaseII) {
                              pdf = (gain == 1) ? pdf_hi : pdf_lo;
                            } else {
                              pdf = pdf_PhaseI;
                            }
                          }
 
                          if (gain == 1) {
                            ped = m_chanPed
                              ? m_tileToolNoiseSample->getPed(drawerIdx, channel, gain,
                                              TileRawChannelUnit::ADCcounts, ctx)
                              : m_ped_HG;
 
                            tFit = m_random->Gaus(0., m_gausC2C);
                            double deformatedTime = m_random->Gaus(m_imperfectionMean, m_imperfectionRms);
                            m_ps[gain]->scalePulse(deformatedTime, deformatedTime);
 
                            // Shift the stored pulses to simulate consecutive BC
                            m_PUAmp[ros-1][drawer][channel].pop_back();
                            m_PUAmp[ros-1][drawer][channel].insert(m_PUAmp[ros-1][drawer][channel].begin(), 0);
 
                            // m_sample_tru is the amplitude for the central BC
                            m_sample_tru = m_PUAmp[ros-1][drawer][channel][(m_nPul - 1) / 2];
 
                            // Fill the new BC at the front
                            addPileUpSample(gain, ros, drawer, channel);
                            m_PUAmp[ros-1][drawer][channel].front() += n_inTimeAmp; // Add amplitude from in-time pulse to true-amp vector
 
                            sample = m_tsg->fillSample(tFit, ped,
                                           m_PUAmp[ros-1][drawer][channel],
                                           pdf, m_gaussNoise,
                                           (int)m_PUAmp[ros-1][drawer][channel].size(),
                                           gain);
                          } else {
                            // Low gain
                            ped = m_chanPed
                              ? m_tileToolNoiseSample->getPed(drawerIdx, channel, gain,
                                              TileRawChannelUnit::ADCcounts, ctx)
                              : m_ped_LG;
 
                            sample = m_tsg->fillSample(tFit, ped,
                                           m_PUAmp[ros-1][drawer][channel],
                                           pdf, m_gaussNoise,
                                           (int)m_PUAmp[ros-1][drawer][channel].size(),
                                           gain);
 
                            // Scale truth amplitude from HG to LG
                            if (m_PhaseII) {
                              m_sample_tru /= 40.0;
                            } else {
                              m_sample_tru /= 64.0;
                            }
                          }
 
                          // Clip the sample
                          if (sample < 0.0) {
                            sample = 0.0;
                          } else if (sample >= m_i_ADCmax) {
                            isHGSaturated = true;
                            if (m_bigain) sample = m_i_ADCmax;
                          }
 
                          samples.clear();
                          samples.push_back(sample);
 
                          ATH_MSG_VERBOSE("New chain ADC " << ros << "/" << drawer << "/" << channel << " - saving gain " << gain);
 
                          TileDigits* digit = tileDigitsPool.nextElementPtr();
                          *digit = TileDigits(m_tileHWID->adc_id(ros, drawer, channel, gain),
                                      std::move(samples));
                          ATH_CHECK(digitsContainer->push_back(digit));
 
                          // Use the “truth” amplitude for rawChannel
                          auto rawChannel = std::make_unique<TileRawChannel>(
                                                     digit->adc_HWID(),
                                                     m_sample_tru,
                                                     tFit,
                                                     m_ootAmp,
                                                     ped);
                          ATH_CHECK(rawChannelContainer->push_back(std::move(rawChannel)));
 
                          // If not bigain, break if HG didn't saturate
                          if (!m_bigain && !isHGSaturated) {
                            break;
                          }
                    }
                }
            }
        }
    }
 
        SG::WriteHandle<TileRawChannelContainer> rawChannelCnt(m_rawChannelContainerKey, ctx);
        ATH_CHECK( rawChannelCnt.record(std::move(rawChannelContainer)) );
 
        SG::WriteHandle<TileDigitsContainer> digitsCnt(m_digitsContainerKey, ctx);
        ATH_CHECK( digitsCnt.record(std::move(digitsContainer)) );
 
    if (!m_simQIE) {
        //delete pdf;
        delete pdf_PhaseI;
        delete pdf_hi;
        delete pdf_lo;
    }
 
    ATH_MSG_DEBUG("Execution completed");
 
    return StatusCode::SUCCESS;
}
 
StatusCode TileDigitsFromPulse::finalize() {
    ATH_MSG_DEBUG("in finalize()");
    delete m_buf;
    delete m_tsg;
    
 
    if (m_useItADist)
        m_itFile->Close();
    if (m_useOotADist)
        m_ootFile->Close();
    if (m_simPUwPoisson)
        m_pileup_AmpDistFile->Close();
 
    return StatusCode::SUCCESS;
}
 
bool TileDigitsFromPulse::makeDist(TFile*& file, TH1F*& hist, const std::string& fileName, const std::string& histName) {
    file = new TFile(fileName.c_str());
    if (file->IsZombie()) {
        ATH_MSG_FATAL("Error reading amplitude distribution from " << fileName << ".");
        return kFALSE;
    }
    TKey *key = file->FindKey(histName.c_str());
    if (key == 0) {
        ATH_MSG_FATAL("Could not find histogram " << histName << " in file " << fileName << ".");
        return kFALSE;
    }
    hist = (TH1F*) file->Get(histName.c_str());
    for (int i = 0; i < m_AmpDistLowLim; i++)
        hist->SetBinContent(i, 0.); // Puts a cut on the amplitude distribution.
    return kTRUE;
 
}
 
bool TileDigitsFromPulse::makeDist(TFile*& file, std::vector<std::vector<TH1F*>>& hists, const std::string& fileName) {
 
    std::string histName;
    TKey *key;
    TH1F* hist;
 
    file = new TFile(fileName.c_str());
    if (file->IsZombie()) {
      ATH_MSG_FATAL("Error reading amplitude distributions from " << fileName << ".");
      return kFALSE;
    }
 
    for(int ros=0; ros<4; ros++){
 
      hists.push_back(std::vector<TH1F*>());
      for(int channel=0; channel<48; channel++){
 
        histName = "ene_ros_" + std::to_string(ros+1) + "_channel_" + std::to_string(channel+1);
 
        key = file->FindKey(histName.c_str());
        if (key == 0) {
          ATH_MSG_FATAL("Could not find histogram " << histName << " in file " << fileName << ".");
          return kFALSE;
        }
 
        hist = (TH1F*) file->Get(histName.c_str());
 
        for (int i = 0; i < m_AmpDistLowLim; i++)
          hist->SetBinContent(i, 0.); // Puts a cut on the amplitude distribution.
 
        hists[ros].push_back(hist);
        hist->Clear();
      }
    }
 
    return kTRUE;
}
 
void TileDigitsFromPulse::addPileUp(double &n_inTimeAmp, int gain, int ros, int drawer, int channel) {
        //Pileup samples
        double amp_1;
        double amp_2;
        double mu = 0; //Interaction per bunch crossing for PU simulation
        if(gain == 1){
                for (int i = 0; i <= m_nPul_eff; i++) {
                        if (((i * 25) % m_BunchSpacing) == 0) {
                                if(m_simPUwPoisson){
                                        mu=m_random->Poisson(m_avgMuForPU);
                                        ATH_MSG_VERBOSE("Effective pulse number " << i);
                                        ATH_MSG_VERBOSE("Number of interactions for simulation: " << mu );
 
                                        for (int imu = 0; imu<mu; imu++){
 
                                                if (m_random->Rndm() < m_pileUpFraction){
                                                        amp_1 = m_pileup_AmpDists[ros-1][channel]->GetRandom();
                                                }
                                                else{
                                                        amp_1 = 0;
                                                }
 
                                                if (m_random->Rndm() < m_pileUpFraction){
                                                        amp_2 = m_pileup_AmpDists[ros-1][channel]->GetRandom();
                                                }
                                                else{
                                                        amp_2 = 0;
                                                }
 
                                                ATH_MSG_VERBOSE("Random amplitudes for PU: " << amp_1 << " " << amp_2);
                                                if(i==0){
                                                        m_PUAmp[ros-1][drawer][channel][m_nPul_eff] += amp_1;
                                                }
                                                else{
                                                        m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] += amp_1;
                                                        m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i] += amp_2;
                                                }
                                        }
 
                                        if(m_PUAmp[ros-1][drawer][channel][m_nPul_eff] < 0) m_PUAmp[ros-1][drawer][channel][m_nPul_eff] = 0;
                                        if(m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] < 0) m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] = 0;
                                        if(m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i] < 0) m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i] = 0;
 
                                        ATH_MSG_VERBOSE("Final amplitudes for pulse " << m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] << " " << m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i]);
                                }
                                else{
                                        m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] = m_useOotADist ? m_ootDist->GetRandom() : m_ootAmp;
                                        m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i] = m_useOotADist ? m_ootDist->GetRandom() : m_ootAmp;
                                }
 
                if(m_simPulseChain){ // Special treatment for pulse-chain simulation, add in-time pulses when initializing true-amp vector
                  if (m_random->Rndm() < m_itAPulseProb){
                    amp_1 = m_useItADist ? m_itDist->GetRandom() : m_inTimeAmp;
                  } else{
                    amp_1 = 0;
                  }
 
                  if (m_random->Rndm() < m_itAPulseProb){
                    amp_2 = m_useItADist ? m_itDist->GetRandom() : m_inTimeAmp;
                  } else{
                    amp_2 = 0;
                  }
 
                  if(i==0){
                    m_PUAmp[ros-1][drawer][channel][m_nPul_eff] += amp_1;
                  }
                  else{
                    m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] += amp_1;
                    m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i] += amp_2;
                  }
                }
                        } else {
                                m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] = 0;
                                m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i] = 0;
                        }
                }
        }
        else if (gain == 0)
        {
                double scaleFactor = m_PhaseII ? 40 : 64;
                n_inTimeAmp /= scaleFactor;
                for (int i = 0; i <= m_nPul_eff; i++) {
                        m_PUAmp[ros-1][drawer][channel][m_nPul_eff + i] /= scaleFactor;
                        m_PUAmp[ros-1][drawer][channel][m_nPul_eff - i] /= scaleFactor;
                }
        }
}
 
void TileDigitsFromPulse::addPileUpSample(int gain, int ros, int drawer, int channel) {
 
    double amp = 0;
    double mu = 0; //Interaction per bunch crossing for PU simulation
 
    auto random = std::make_unique<TRandom3>(m_seed);
 
    mu=random->Poisson(m_avgMuForPU);
 
    if(gain == 1){
        for (int imu = 0; imu<mu; imu++){
            if (random->Rndm() < m_pileUpFraction){
                amp = m_pileup_AmpDists[ros-1][channel]->GetRandom();
            }
            else{
                amp = 0;
            }
 
            m_PUAmp[ros-1][drawer][channel].front() += amp;
        }
    }
    else if (gain == 0)
    {
        if(!m_PhaseII){
            m_PUAmp[ros-1][drawer][channel].front() /= 64;
        }
        else{
            m_PUAmp[ros-1][drawer][channel].front() /= 40;
        }
    }
 
    if(m_PUAmp[ros-1][drawer][channel].front() < 0) m_PUAmp[ros-1][drawer][channel].front() = 0;
}