diff --git a/Detectors/FIT/FT0/base/include/FT0Base/FT0DigParam.h b/Detectors/FIT/FT0/base/include/FT0Base/FT0DigParam.h index 2bf774859aa22..074d91bb04b27 100644 --- a/Detectors/FIT/FT0/base/include/FT0Base/FT0DigParam.h +++ b/Detectors/FIT/FT0/base/include/FT0Base/FT0DigParam.h @@ -31,8 +31,8 @@ struct FT0DigParam : o2::conf::ConfigurableParamHelper { 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 = 40; // Tclu units (40 for pp and 1433 for PbPb in Run3) + int mtrg_semicentral_trh = 20; // Tclu units (20 for pp and 35 for PbPb in Run3) float mMip_in_V = 7; // MIP to mV float mPe_in_mip = 0.004; // invserse Np.e. in MIP 1./250. @@ -43,11 +43,12 @@ struct FT0DigParam : o2::conf::ConfigurableParamHelper { 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 - float mMV_2_Nchannels = 2.2857143; // amplitude channel 7 mV ->16channels - float mMV_2_NchannelsInverse = 0.437499997; // inverse amplitude channel 7 mV ->16channels + float mMV_2_Nchannels = 2.; // amplitude channel 7 mV ->14channels + float mMV_2_NchannelsInverse = 0.5; // inverse amplitude channel 7 mV ->14channels (nowhere used) O2ParamDef(FT0DigParam, "FT0DigParam"); }; diff --git a/Detectors/FIT/FT0/simulation/src/Digitizer.cxx b/Detectors/FIT/FT0/simulation/src/Digitizer.cxx index aca012f1bc5a9..de432a85765c7 100644 --- a/Detectors/FIT/FT0/simulation/src/Digitizer.cxx +++ b/Detectors/FIT/FT0/simulation/src/Digitizer.cxx @@ -16,6 +16,13 @@ #include "CommonConstants/PhysicsConstants.h" #include "CommonDataFormat/InteractionRecord.h" +#include "DataFormatsFT0/LookUpTable.h" +#include "FT0Base/Constants.h" +#include +#include +#include +#include + #include "TMath.h" #include "TRandom.h" #include @@ -35,24 +42,84 @@ namespace o2::ft0 template 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 +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) { @@ -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 mChID2PMhash{}; + static std::map mMapPMhash2isAside; // hashed PM -> is A side + + if (!pmLutInitialized) { + std::map 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 sum_ampl_ipmt(nPMTs, 0); + // Per-PM summed charge (like in digits2trgFT0) + std::map 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; @@ -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(ipmt)]] += static_cast(amp); + } } digitsCh.emplace_back(ipmt, smeared_time, int(amp), chain); nStored++; @@ -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; @@ -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; diff --git a/Detectors/FIT/FV0/simulation/include/FV0Simulation/FV0DigParam.h b/Detectors/FIT/FV0/simulation/include/FV0Simulation/FV0DigParam.h index 383fa4cb494c1..6462323a279b7 100644 --- a/Detectors/FIT/FV0/simulation/include/FV0Simulation/FV0DigParam.h +++ b/Detectors/FIT/FV0/simulation/include/FV0Simulation/FV0DigParam.h @@ -69,6 +69,10 @@ struct FV0DigParam : o2::conf::ConfigurableParamHelper { 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 = 2; // trigger Nchannels + float InnerChargeLevel = 4; // InnerRingsChargeLevel + float OuterChargeLevel = 4; // OuterRingsChargeLevel + float ChargeLevel = 8; // ChargeLevel O2ParamDef(FV0DigParam, "FV0DigParam"); }; diff --git a/Detectors/FIT/FV0/simulation/src/Digitizer.cxx b/Detectors/FIT/FV0/simulation/src/Digitizer.cxx index 8c1d2dc8824e2..3237f9bab7879 100644 --- a/Detectors/FIT/FV0/simulation/src/Digitizer.cxx +++ b/Detectors/FIT/FV0/simulation/src/Digitizer.cxx @@ -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++) { @@ -149,8 +149,8 @@ void Digitizer::process(const std::vector& 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, @@ -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]); @@ -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; @@ -285,8 +287,10 @@ void Digitizer::storeBC(const BCCache& bc, avgTime += iCfdZero; if (iPmt < 24) { nSignalInner++; + totalChargeInnerRing += iTotalCharge; } else { nSignalOuter++; + totalChargeOuterRing += iTotalCharge; } } } @@ -300,13 +304,15 @@ 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 + isAIn = 0.125 * totalChargeInnerRing > 2 * FV0DigParam::Instance().InnerChargeLevel; // ring 1,2 and 3 + // isAOut = nSignalOuter > FV0DigParam::Instance().NchannelsLevel; // ring 4 and 5 + isAOut = 0.125 * totalChargeOuterRing > 2 * FV0DigParam::Instance().OuterChargeLevel; // ring 4 and 5 + isTotalCharge = 0.125 * totalChargeAllRing > 2 * FV0DigParam::Instance().ChargeLevel; Triggers triggers; const int unusedCharge = o2::fit::Triggers::DEFAULT_AMP; @@ -314,10 +320,10 @@ void Digitizer::storeBC(const BCCache& bc, 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); } @@ -342,8 +348,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; } @@ -400,7 +406,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); }