Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions Detectors/FIT/FT0/base/include/FT0Base/FT0DigParam.h
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,8 @@ struct FT0DigParam : o2::conf::ConfigurableParamHelper<FT0DigParam> {
float mAmpRecordUp = 15; // to [ns]
float hitTimeOffsetA = 0; ///< hit time offset on the A side [ns]
float hitTimeOffsetC = 0; ///< hit time offset on the C side [ns]
int mtrg_central_trh = 600.; // channels
int mtrg_semicentral_trh = 300.; // channels
int mtrg_central_trh = 1433.; // Tclu
int mtrg_semicentral_trh = 35.; // Tclu units = TCM charge level unit = 16 ADCunits

float mMip_in_V = 7; // MIP to mV
float mPe_in_mip = 0.004; // invserse Np.e. in MIP 1./250.
Expand All @@ -43,6 +43,7 @@ struct FT0DigParam : o2::conf::ConfigurableParamHelper<FT0DigParam> {
float mNoiseVar = 0.1; // noise level
float mNoisePeriod = 1 / 0.9; // GHz low frequency noise period;
short mTime_trg_gate = 153; // #channels as in TCM as in Pilot beams ('OR gate' setting in TCM tab in ControlServer)
short mTime_trg_vertex_gate = 100; // #channels as in TCM as in Pilot beams ('OR gate' setting in TCM tab in ControlServer)
float mAmpThresholdForReco = 5; // only channels with amplitude higher will participate in calibration and collision time: 0.3 MIP
short mTimeThresholdForReco = 1000; // only channels with time below will participate in calibration and collision time

Expand Down
185 changes: 172 additions & 13 deletions Detectors/FIT/FT0/simulation/src/Digitizer.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,13 @@
#include "CommonConstants/PhysicsConstants.h"
#include "CommonDataFormat/InteractionRecord.h"

#include "DataFormatsFT0/LookUpTable.h"
#include "FT0Base/Constants.h"
#include <map>
#include <array>
#include <regex>
#include <string>

#include "TMath.h"
#include "TRandom.h"
#include <TH1F.h>
Expand All @@ -35,24 +42,84 @@ namespace o2::ft0
template <typename Float>
Float signalForm_i(Float x)
{
using namespace std;
Float const a = -0.45458;
Float const b = -0.83344945;
return x > Float(0) ? -(exp(b * x) - exp(a * x)) / Float(7.8446501) : Float(0);
float p0, p1, p2, p3, p4, p5, p6, p7;
p0 = 1.30853;
p1 = 0.516807;
p2 = 3.36714;
p3 = -1.01206;
p4 = 1.42832;
p5 = 1.1589;
p6 = 1.22019;
p7 = 0.426818;

Double_t val = 0;

if (x > p3) {
Double_t x1 = x - p3;
Double_t arg1 = (log(x1) - p0) / p1;
val += p2 * (1.0 / (x1 * p1 * sqrt(2 * TMath::Pi()))) * exp(-0.5 * arg1 * arg1);
}

if (x > p7) {
Double_t x2 = x - p7;
Double_t arg2 = (log(x2) - p4) / p5;
val += p6 * (1.0 / (x2 * p5 * sqrt(2 * TMath::Pi()))) * exp(-0.5 * arg2 * arg2);
}

return val;
};

// integrated signal shape function
inline float signalForm_integral(float x)
{
using namespace std;
double const a = -0.45458;
double const b = -0.83344945;
if (x < 0) {
x = 0;
float p0, p1, p2, p3, p4, p5, p6, p7;
p0 = 1.30853;
p1 = 0.516807;
p2 = 3.36714;
p3 = -1.01206;
p4 = 1.42832;
p5 = 1.1589;
p6 = 1.22019;
p7 = 0.426818;
Double_t val = 0;

if (x > p3) {
Double_t x1 = x - p3;
Double_t z1 = (log(x1) - p0) / (sqrt(2) * p1);
val += p2 * 0.5 * (1 + TMath::Erf(z1)); // norm1 * CDF1
}
return -(exp(b * x) / b - exp(a * x) / a) / 7.8446501;

if (x > p7) {
Double_t x2 = x - p7;
Double_t z2 = (log(x2) - p4) / (sqrt(2) * p5);
val += p6 * 0.5 * (1 + TMath::Erf(z2)); // norm2 * CDF2
}

return val;
};
/*
// signal shape function
template <typename Float>
Float signalForm_i(Float x)
{
using namespace std;
Float const a = -0.45458;
Float const b = -0.83344945;
return x > Float(0) ? -(exp(b * x) - exp(a * x)) / Float(7.8446501) : Float(0);
};

// integrated signal shape function
inline float signalForm_integral(float x)
{
using namespace std;
double const a = -0.45458;
double const b = -0.83344945;
if (x < 0) {
x = 0;
}
return -(exp(b * x) / b - exp(a * x) / a) / 7.8446501;
};
*/
// SIMD version of the integrated signal shape function
inline Vc::float_v signalForm_integralVc(Vc::float_v x)
{
Expand Down Expand Up @@ -249,8 +316,64 @@ void Digitizer::storeBC(BCCache& bc,
if (bc.hits.empty()) {
return;
}
// Initialize mapping channelID -> PM hash and PM side (A/C) using FT0 LUT
static bool pmLutInitialized = false;
static std::array<uint8_t, o2::ft0::Constants::sNCHANNELS_PM> mChID2PMhash{};
static std::map<uint8_t, bool> mMapPMhash2isAside; // hashed PM -> is A side

if (!pmLutInitialized) {
std::map<std::string, uint8_t> mapFEE2hash; // module name -> hashed PM id
uint8_t tcmHash = 0;

const auto& lut = o2::ft0::SingleLUT::Instance().getVecMetadataFEE();
auto lutSorted = lut;
std::sort(lutSorted.begin(), lutSorted.end(),
[](const auto& first, const auto& second) { return first.mModuleName < second.mModuleName; });

uint8_t binPos = 0;
for (const auto& lutEntry : lutSorted) {
const auto& moduleName = lutEntry.mModuleName;
const auto& moduleType = lutEntry.mModuleType;
const auto& strChID = lutEntry.mChannelID;

auto [it, inserted] = mapFEE2hash.insert({moduleName, binPos});
if (inserted) {
if (moduleName.find("PMA") != std::string::npos) {
mMapPMhash2isAside.insert({binPos, true});
} else if (moduleName.find("PMC") != std::string::npos) {
mMapPMhash2isAside.insert({binPos, false});
}
++binPos;
}

if (std::regex_match(strChID, std::regex("^[0-9]{1,3}$"))) {
int chID = std::stoi(strChID);
if (chID < o2::ft0::Constants::sNCHANNELS_PM) {
mChID2PMhash[chID] = mapFEE2hash[moduleName];
} else {
LOG(fatal) << "Incorrect LUT entry: chID " << strChID << " | " << moduleName;
}
} else if (moduleType != "TCM") {
LOG(fatal) << "Non-TCM module w/o numerical chID: chID " << strChID << " | " << moduleName;
} else { // TCM
tcmHash = mapFEE2hash[moduleName];
}
}

pmLutInitialized = true;
}

int n_hit_A = 0, n_hit_C = 0, mean_time_A = 0, mean_time_C = 0;
int summ_ampl_A = 0, summ_ampl_C = 0;
int sum_A_ampl = 0, sum_C_ampl = 0;
int nPMTs = mGeometry.NCellsA * 4 + mGeometry.NCellsC * 4;
std::vector<int> sum_ampl_ipmt(nPMTs, 0);
// Per-PM summed charge (like in digits2trgFT0)
std::map<uint8_t, int> mapPMhash2sumAmpl;
for (const auto& entry : mMapPMhash2isAside) {
mapPMhash2sumAmpl.insert({entry.first, 0});
}

int vertex_time;
const auto& params = FT0DigParam::Instance();
int first = digitsCh.size(), nStored = 0;
Expand Down Expand Up @@ -297,6 +420,10 @@ void Digitizer::storeBC(BCCache& bc,
if (is_time_in_signal_gate) {
chain |= (1 << o2::ft0::ChannelData::EEventDataBit::kIsCFDinADCgate);
chain |= (1 << o2::ft0::ChannelData::EEventDataBit::kIsEventInTVDC);
// Sum channel charge per PM (similar logic as in digits2trgFT0)
if (ipmt < o2::ft0::Constants::sNCHANNELS_PM) {
mapPMhash2sumAmpl[mChID2PMhash[static_cast<uint8_t>(ipmt)]] += static_cast<int>(amp);
}
}
digitsCh.emplace_back(ipmt, smeared_time, int(amp), chain);
nStored++;
Expand All @@ -308,6 +435,8 @@ void Digitizer::storeBC(BCCache& bc,
continue;
}

sum_ampl_ipmt[ipmt] += amp;

if (is_A_side) {
n_hit_A++;
summ_ampl_A += amp;
Expand All @@ -318,17 +447,47 @@ void Digitizer::storeBC(BCCache& bc,
mean_time_C += smeared_time;
}
}

for (size_t i = 0; i < sum_ampl_ipmt.size(); i++) {
sum_ampl_ipmt[i] = sum_ampl_ipmt[i] >> 3;
if (i < 4 * mGeometry.NCellsA) {
sum_A_ampl += sum_ampl_ipmt[i];
} else {
sum_C_ampl += sum_ampl_ipmt[i];
}
}

// Sum over PMs (using per-PM map) for debug/monitoring
int sum_PM_ampl_debug = 0;
int sum_PM_ampl_A_debug = 0;
int sum_PM_ampl_C_debug = 0;
for (const auto& entry : mapPMhash2sumAmpl) {
int pmAmpl = (entry.second >> 3);
sum_PM_ampl_debug += pmAmpl;
auto itSide = mMapPMhash2isAside.find(entry.first);
if (itSide != mMapPMhash2isAside.end()) {
if (itSide->second) {
sum_PM_ampl_A_debug += pmAmpl;
} else {
sum_PM_ampl_C_debug += pmAmpl;
}
}
}
LOG(debug) << "Sum PM amplitude (LUT-based): total=" << sum_PM_ampl_debug
<< " A-side=" << sum_PM_ampl_A_debug
<< " C-side=" << sum_PM_ampl_C_debug;

Bool_t is_A, is_C, isVertex, is_Central, is_SemiCentral = 0;
is_A = n_hit_A > 0;
is_C = n_hit_C > 0;
is_Central = summ_ampl_A + summ_ampl_C >= params.mtrg_central_trh;
is_SemiCentral = summ_ampl_A + summ_ampl_C >= params.mtrg_semicentral_trh;
is_Central = sum_PM_ampl_A_debug + sum_PM_ampl_C_debug >= 2 * params.mtrg_central_trh;
is_SemiCentral = sum_PM_ampl_A_debug + sum_PM_ampl_C_debug >= 2 * params.mtrg_semicentral_trh && !is_Central;
uint32_t amplA = is_A ? summ_ampl_A * 0.125 : -5000; // sum amplitude A side / 8 (hardware)
uint32_t amplC = is_C ? summ_ampl_C * 0.125 : -5000; // sum amplitude C side / 8 (hardware)
int timeA = is_A ? mean_time_A / n_hit_A : -5000; // average time A side
int timeC = is_C ? mean_time_C / n_hit_C : -5000; // average time C side
vertex_time = (timeC - timeA) * 0.5;
isVertex = is_A && is_C && (vertex_time > -params.mTime_trg_gate && vertex_time < params.mTime_trg_gate);
isVertex = is_A && is_C && (vertex_time > -params.mTime_trg_vertex_gate && vertex_time < params.mTime_trg_vertex_gate);
LOG(debug) << " A " << is_A << " timeA " << timeA << " mean_time_A " << mean_time_A << " n_hit_A " << n_hit_A << " C " << is_C << " timeC " << timeC << " mean_time_C " << mean_time_C << " n_hit_C " << n_hit_C << " vertex_time " << vertex_time;
Triggers triggers;
bool isLaser = false;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,10 @@ struct FV0DigParam : o2::conf::ConfigurableParamHelper<FV0DigParam> {
uint8_t defaultChainQtc = 0x48; // only 2 flags are set by default in simulation: kIsCFDinADCgate and kIsEventInTVDC
float mAmpThresholdForReco = 24; // only channels with amplitude higher will participate in calibration and collision time
short mTimeThresholdForReco = 1000; // only channels with time below will participate in calibration and collision time
int NchannelsLevel = 8; // trigger Nchannels paramter
float InnerRingsChargeLevel = 0; // InnerRingsChargeLevel, Tclu
float OuterRingsChargeLevel = 0; // TcluOuterRingsChargeLevel, Tclu
float ChargeLevel = 10; // ChargeLevel, Tclu

O2ParamDef(FV0DigParam, "FV0DigParam");
};
Expand Down
36 changes: 20 additions & 16 deletions Detectors/FIT/FV0/simulation/src/Digitizer.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,8 @@ void Digitizer::clear()
void Digitizer::init()
{
LOG(info) << "init";
mNBins = FV0DigParam::Instance().waveformNbins; //Will be computed using detector set-up from CDB
mBinSize = FV0DigParam::Instance().waveformBinWidth; //Will be set-up from CDB
mNBins = FV0DigParam::Instance().waveformNbins; // Will be computed using detector set-up from CDB
mBinSize = FV0DigParam::Instance().waveformBinWidth; // Will be set-up from CDB
mNTimeBinsPerBC = std::lround(o2::constants::lhc::LHCBunchSpacingNS / mBinSize); // 1920 bins/BC

for (Int_t detID = 0; detID < Constants::nFv0Channels; detID++) {
Expand Down Expand Up @@ -149,8 +149,8 @@ void Digitizer::process(const std::vector<o2::fv0::Hit>& hits,

createPulse(mipFraction, hit.GetTrackID(), hitTime, hit.GetPos().R(), cachedIR, nCachedIR, detId);

} //while loop
} //hitloop
} // while loop
} // hitloop
}

void Digitizer::createPulse(float mipFraction, int parID, const double hitTime, const float hitR,
Expand Down Expand Up @@ -200,7 +200,7 @@ void Digitizer::createPulse(float mipFraction, int parID, const double hitTime,
}
added[ir] = true;
}
///Add MC labels to BCs for those contributed to the PMT signal
/// Add MC labels to BCs for those contributed to the PMT signal
for (int ir = 0; ir < nCachedIR; ir++) {
if (added[ir]) {
auto bcCache = getBCCache(cachedIR[ir]);
Expand Down Expand Up @@ -238,6 +238,8 @@ void Digitizer::storeBC(const BCCache& bc,
int8_t nTotFiredCells = 0;
int8_t nTrgFiredCells = 0; // number of fired cells, that follow additional trigger conditions (time gate)
int totalChargeAllRing = 0;
int totalChargeInnerRing = 0;
int totalChargeOuterRing = 0;
int32_t avgTime = 0;
double nSignalInner = 0;
double nSignalOuter = 0;
Expand Down Expand Up @@ -285,8 +287,10 @@ void Digitizer::storeBC(const BCCache& bc,
avgTime += iCfdZero;
if (iPmt < 24) {
nSignalInner++;
totalChargeInnerRing += iTotalCharge;
} else {
nSignalOuter++;
totalChargeOuterRing += iTotalCharge;
}
}
}
Expand All @@ -300,24 +304,24 @@ void Digitizer::storeBC(const BCCache& bc,
} else {
avgTime = o2::fit::Triggers::DEFAULT_TIME;
}
///Triggers for FV0
bool isA, isAIn, isAOut, isCen, isSCen;
/// Triggers for FV0
bool isA, isNchannels, isAIn, isAOut, isTotalCharge;
isA = nTrgFiredCells > 0;
isAIn = nSignalInner > 0; // ring 1,2 and 3
isAOut = nSignalOuter > 0; // ring 4 and 5
isCen = totalChargeAllRing > FV0DigParam::Instance().adcChargeCenThr;
isSCen = totalChargeAllRing > FV0DigParam::Instance().adcChargeSCenThr;
isNchannels = nTrgFiredCells > FV0DigParam::Instance().NchannelsLevel;
isAIn = nSignalInner > FV0DigParam::Instance().NchannelsLevel; // ring 1,2 and 3
isAOut = nSignalOuter > FV0DigParam::Instance().NchannelsLevel; // ring 4 and 5
isTotalCharge = 0.125 * totalChargeAllRing > 2 * FV0DigParam::Instance().ChargeLevel;

Triggers triggers;
const int unusedCharge = o2::fit::Triggers::DEFAULT_AMP;
const int unusedTime = o2::fit::Triggers::DEFAULT_TIME;
const int unusedZero = o2::fit::Triggers::DEFAULT_ZERO;
const bool unusedBitsInSim = false; // bits related to laser and data validity
const bool bitDataIsValid = true;
triggers.setTriggers(isA, isAIn, isAOut, isCen, isSCen, nTrgFiredCells, (int8_t)unusedZero,
triggers.setTriggers(isA, isAIn, isAOut, isTotalCharge, isNchannels, nTrgFiredCells, (int8_t)unusedZero,
(int32_t)(0.125 * totalChargeAllRing), (int32_t)unusedCharge, (int16_t)avgTime, (int16_t)unusedTime, unusedBitsInSim, unusedBitsInSim, bitDataIsValid);
digitsBC.emplace_back(first, nTotFiredCells, bc, triggers, mEventId - 1);
digitsTrig.emplace_back(bc, isA, isAIn, isAOut, isCen, isSCen);
digitsTrig.emplace_back(bc, isA, isAIn, isAOut, isTotalCharge, isNchannels);
for (auto const& lbl : bc.labels) {
labels.addElement(nBC, lbl);
}
Expand All @@ -342,8 +346,8 @@ Int_t Digitizer::SimulateLightYield(Int_t pmt, Int_t nPhot) const
//---------------------------------------------------------------------------
Float_t Digitizer::IntegrateCharge(const ChannelDigitF& pulse) const
{
int const chargeIntMin = FV0DigParam::Instance().isIntegrateFull ? 0 : (FV0DigParam::Instance().avgCfdTimeForMip - 6.0) / mBinSize; //Charge integration offset (cfd mean time - 6 ns)
int const chargeIntMax = FV0DigParam::Instance().isIntegrateFull ? mNTimeBinsPerBC : (FV0DigParam::Instance().avgCfdTimeForMip + 14.0) / mBinSize; //Charge integration offset (cfd mean time + 14 ns)
int const chargeIntMin = FV0DigParam::Instance().isIntegrateFull ? 0 : (FV0DigParam::Instance().avgCfdTimeForMip - 6.0) / mBinSize; // Charge integration offset (cfd mean time - 6 ns)
int const chargeIntMax = FV0DigParam::Instance().isIntegrateFull ? mNTimeBinsPerBC : (FV0DigParam::Instance().avgCfdTimeForMip + 14.0) / mBinSize; // Charge integration offset (cfd mean time + 14 ns)
if (chargeIntMin < 0 || chargeIntMin > mNTimeBinsPerBC || chargeIntMax > mNTimeBinsPerBC) {
LOG(fatal) << "invalid indicess: chargeInMin=" << chargeIntMin << " chargeIntMax=" << chargeIntMax;
}
Expand Down Expand Up @@ -400,7 +404,7 @@ float Digitizer::getDistFromCellCenter(UInt_t cellId, double hitx, double hity)
double a = -(y0 - pCell->y) / (x0 - pCell->x);
double b = 1;
double c = -(y0 - a * x0);
//Return the distance from hit to this line
// Return the distance from hit to this line
return (a * hitx + b * hity + c) / TMath::Sqrt(a * a + b * b);
}

Expand Down
Loading