BmnGemStripModule.cxx

Go to the documentation of this file.

#include "BmnGemStripModule.h"
 
#include "BmnGemStripDigitizer.h"   //AZ-110124
#include "CbmMCTrack.h"             //AZ-110124
#include "CbmStsPoint.h"            //AZ-110124
#include "FairLogger.h"
#include "FairRunSim.h"   //AZ-110124
#include "TF1.h"
#include "TH1F.h"
#include "TRandom.h"
 
#include <TDatabasePDG.h>   //AZ-110124
#include <TParticlePDG.h>   //AZ-110124
#include <algorithm>
#include <atomic>               // [PARALLEL-STD]
#include <condition_variable>   // [PARALLEL-STD]
#include <functional>           // [PARALLEL-STD]
#include <limits>               // [PARALLEL-STD]
#include <map>
#include <mutex>
#include <random>   // [PARALLEL-STD]
#include <thread>   // [PARALLEL-STD]
 
// ANN models: own inference (MLP) and frugally-deep (fdeep is used only in BmnGemStripFdeepModel.cxx)
#include "BmnGemStripFdeepModel.h"
#include "BmnGemStripMLPModel.h"
//------------------------------------------------------------------------------
 
// ANN generation of the clusters: static bool for all modules
Bool_t BmnGemStripModule::fUseANNGeneration = kFALSE;
 
// Path to an ANN model
void BmnGemStripModule::SetANNModelPath(TString path)
{
    BmnGemStripMLPModel::SetDefaultPath(path.Data());   // "" -> the default path of the MLP model
}
 
TString BmnGemStripModule::GetANNModelPath()
{
    return TString(BmnGemStripMLPModel::GetDefaultPath().c_str());
}
 
BmnGemStripModule::BmnGemStripModule()
{
    Verbosity = false;
 
    ZStartModulePosition = 0.0;
 
    DriftGapThickness = 0.3;            // cm
    FirstTransferGapThickness = 0.25;   // cm
    SecondTransferGapThickness = 0.2;   // cm
    InductionGapThickness = 0.15;       // cm
    ModuleThickness =
        DriftGapThickness + FirstTransferGapThickness + SecondTransferGapThickness + InductionGapThickness;   // cm
 
    ElectronDriftDirection = ForwardZAxisEDrift;
 
    XMinModule = 0.0;
    XMaxModule = 0.0;
    YMinModule = 0.0;
    YMaxModule = 0.0;
 
    AvalancheGenerationSeed = 0;
    MaxClustersPerCm = 0.0;       // no limit
    FillGenerationData = false;   // no generation data in the strip methods (see SetFillGenerationData)
    AvalancheRadius = 0.10;       // cm
    Gain = 1.0;   // gain level (for each electron signal - in RealPointFull or for strip signal - in RealPointFullOne)
}
 
BmnGemStripModule::BmnGemStripModule(Double_t z_start_pos,
                                     ElectronDriftDirectionInModule edrift_direction,
                                     Double_t DriftGap,
                                     Double_t FTransferGap,
                                     Double_t STransferGap,
                                     Double_t InductionGap)
{
    Verbosity = false;
 
    ZStartModulePosition = z_start_pos;
 
    DriftGapThickness = DriftGap;                // cm
    FirstTransferGapThickness = FTransferGap;    // cm
    SecondTransferGapThickness = STransferGap;   // cm
    InductionGapThickness = InductionGap;        // cm
    ModuleThickness =
        DriftGapThickness + FirstTransferGapThickness + SecondTransferGapThickness + InductionGapThickness;   // cm
 
    ElectronDriftDirection = edrift_direction;
 
    XMinModule = 0.0;
    XMaxModule = 0.0;
    YMinModule = 0.0;
    YMaxModule = 0.0;
 
    AvalancheGenerationSeed = 0;
    MaxClustersPerCm = 0.0;       // no limit
    FillGenerationData = false;   // no generation data in the strip methods (see SetFillGenerationData)
    AvalancheRadius = 0.10;       // cm
    Gain = 1.0;
}
 
BmnGemStripModule::~BmnGemStripModule() {}
 
//------------------------------------------------------------------------------
 
Double_t BmnGemStripModule::GetZPositionRegistered()
{
    if (ElectronDriftDirection == ForwardZAxisEDrift) {
        return ZStartModulePosition + DriftGapThickness * 0.5;
    } else {
        return ZStartModulePosition + ModuleThickness - DriftGapThickness * 0.5;
    }
}
 
void BmnGemStripModule::AddStripLayer(BmnGemStripLayer strip_layer)
{
    StripLayers.push_back(strip_layer);
 
    DefineModuleBorders();
}
 
Bool_t BmnGemStripModule::SetStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
{
    if (layer_num >= 0 && layer_num < (int)StripLayers.size()) {
        return StripLayers[layer_num].SetStripSignal(strip_num, signal);
    }
    return false;
}
 
Bool_t BmnGemStripModule::AddStripSignalInLayer(Int_t layer_num, Int_t strip_num, Double_t signal)
{
    if (layer_num >= 0 && layer_num < (int)StripLayers.size()) {
        return StripLayers[layer_num].AddStripSignal(strip_num, signal);
    }
    return false;
}
 
Bool_t BmnGemStripModule::SetStripMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch mc_match)
{
    if (layer_num >= 0 && layer_num < (int)StripLayers.size()) {
        return StripLayers[layer_num].SetStripMatch(strip_num, mc_match);
    }
    return false;
}
 
Bool_t BmnGemStripModule::SetStripDigitNumberMatchInLayer(Int_t layer_num, Int_t strip_num, BmnMatch digit_num_match)
{
    if (layer_num >= 0 && layer_num < (int)StripLayers.size()) {
        return StripLayers[layer_num].SetStripDigitNumberMatch(strip_num, digit_num_match);
    }
    return false;
}
 
Double_t BmnGemStripModule::GetStripSignalInLayer(Int_t layer_num, Int_t strip_num)
{
    if (layer_num >= 0 && layer_num < (int)StripLayers.size()) {
        return StripLayers[layer_num].GetStripSignal(strip_num);
    }
    return -1;
}
 
BmnMatch BmnGemStripModule::GetStripMatchInLayer(Int_t layer_num, Int_t strip_num)
{
    if (layer_num >= 0 && layer_num < (int)StripLayers.size()) {
        return StripLayers[layer_num].GetStripMatch(strip_num);
    }
    return BmnMatch();
}
 
BmnMatch BmnGemStripModule::GetStripDigitNumberMatchInLayer(Int_t layer_num, Int_t strip_num)
{
    if (layer_num >= 0 && layer_num < (int)StripLayers.size()) {
        return StripLayers[layer_num].GetStripDigitNumberMatch(strip_num);
    }
    return BmnMatch();
}
 
void BmnGemStripModule::ResetModuleData()
{
    for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
        StripLayers[ilayer].ResetLayer();
    }
 
    DefineModuleBorders();
 
    ResetRealPoints();
    ResetIntersectionPoints();
    ResetPseudoIntersections();
 
    ResetElectronPointsAndClusters();   // test
}
 
Bool_t BmnGemStripModule::IsPointInsideModule(Double_t x, Double_t y)
{
    for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
        if (StripLayers[ilayer].IsPointInsideStripLayer(x, y))
            return true;
    }
    return false;
}
 
Bool_t BmnGemStripModule::IsPointInsideModule(Double_t x, Double_t y, Double_t z)
{
    // old condition: if( z >= ZStartModulePosition && z <= (ZStartModulePosition+ModuleThickness) )
    Double_t coord_eps = 0.01;   // 100 um
    if (((z > ZStartModulePosition) || (fabs(z - ZStartModulePosition) < coord_eps))
        && ((z < (ZStartModulePosition + ModuleThickness))
            || (fabs(z - (ZStartModulePosition + ModuleThickness))) < coord_eps))
    {
        if (IsPointInsideModule(x, y))
            return true;
    }
    return false;
}
 
Bool_t BmnGemStripModule::IsPointInsideZThickness(Double_t z)
{
    Double_t coord_eps = 0.01;   // 100 um
    if (((z > ZStartModulePosition) || (fabs(z - ZStartModulePosition) < coord_eps))
        && ((z < (ZStartModulePosition + ModuleThickness))
            || (fabs(z - (ZStartModulePosition + ModuleThickness))) < coord_eps))
    {
        return true;
    }
    return false;
}
 
// Window size of a strip ANN model: odd (a center strip) and >= 3. Empty string - OK, otherwise - the reason
static std::string CheckANNStripWindowSize(Int_t n_out)
{
    if (n_out < 3 || n_out % 2 == 0)
        return "model output size " + std::to_string(n_out) + ": the strip window must be odd and >= 3";
    return "";
}
 
// Center point of the ANN strip window: the point of the track segment (entry (x, y, z) -> exit (x_out, y_out, z_out))
// at the depth of the middle of the drift gap, where most of the signal comes from (the electrons of the drift gap
// pass all three GEM cascades). The drift gap is at the cathode side: at the module start for the forward drift
// and at the module end for the backward one. The Lorentz shift is not included (the ANN learns it from bfield).
// Tracks that do not reach this depth (they start/stop inside the module or leave it through a side): the nearest
// end of the segment. The same point for pz > 0 and pz < 0 (only the segment matters).
void BmnGemStripModule::GetANNWindowCenterPoint(Double_t x,
                                                Double_t y,
                                                Double_t z,
                                                Double_t x_out,
                                                Double_t y_out,
                                                Double_t z_out,
                                                Double_t& x_c,
                                                Double_t& y_c) const
{
    // distance from the middle of the drift gap to the readout plane
    const Double_t zdist_mid =
        InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness + 0.5 * DriftGapThickness;
    const Double_t z_mid = (ElectronDriftDirection == ForwardZAxisEDrift)
                               ? ZStartModulePosition + ModuleThickness - zdist_mid   // readout at the module end
                               : ZStartModulePosition + zdist_mid;                    // readout at the module start
    Double_t t = (TMath::Abs(z_out - z) > 1e-9) ? (z_mid - z) / (z_out - z) : 0.0;
    t = TMath::Max(0.0, TMath::Min(1.0, t));
    x_c = x + t * (x_out - x);
    y_c = y + t * (y_out - y);
}
//------------------------------------------------------------------------------
 
// Universal method: the standard method or the ANN one (see SetUseANNGeneration())
Bool_t BmnGemStripModule::AddRealPoint(Double_t x,
                                       Double_t y,
                                       Double_t z,
                                       Double_t px,
                                       Double_t py,
                                       Double_t pz,
                                       Double_t signal,
                                       Int_t refID)
{
    if (!fUseANNGeneration)
        return AddRealPointFull(x, y, z, px, py, pz, signal, refID);
    else
        return AddRealPointFull_ANNStripGenerator_MLP_FastStand(x, y, z, px, py, pz, signal, refID, 0, -1);
}
//------------------------------------------------------------------------------
 
// Add point with full realistic behaviour (tracking through the module, avalanche effects)
Bool_t BmnGemStripModule::AddRealPointFull(Double_t x,
                                           Double_t y,
                                           Double_t z,
                                           Double_t px,
                                           Double_t py,
                                           Double_t pz,
                                           Double_t signal,
                                           Int_t refID)
{
 
    // Condition: a start point is inside the module and outside its dead zone (in strip layers)
    if (IsPointInsideModule(x, y)) {
 
        // AZ-260322 gRandom->SetSeed(AvalancheGenerationSeed);
 
        // AZ-110124 - simulation of ionization clusters
        BmnGemStripDigitizer* digi = (BmnGemStripDigitizer*)FairRunSim::Instance()->GetTask("GemDigitizer");
        TClonesArray* points = digi->GetMCPoints();
        TClonesArray* mctracks = digi->GetMCTracks();
        CbmStsPoint* point = (CbmStsPoint*)points->UncheckedAt(refID);
        TVector3 mom;
        point->Momentum(mom);
        CbmMCTrack* mcTr = (CbmMCTrack*)mctracks->UncheckedAt(point->GetTrackID());
        TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(mcTr->GetPdgCode());
        Double_t mnoc = 1.0, mass = 3.7283999;   // He4
        if (part)
            mass = part->Mass();
        if (mass > 0.0001) {
            TLorentzVector lorv;
            lorv.SetVectM(mom, mass);
            Double_t beta = lorv.Beta();
            Double_t gamma = lorv.Gamma();
            Double_t charge = (part) ? part->Charge() / 3 : 1;
            mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
            if (mnoc < 0.5)
                return kFALSE;
            // cout << " ***clus*** " << mass << " " << beta << " " << gamma << " " << mnoc << endl;
 
            // test couts
            if (0) {
                cout << "particle,  pdg_code = " << part->PdgCode() << "\n";
                cout << "  px:py:pz = ( " << px << " : " << py << " : " << pz << " )\n";
                cout << "  mass = " << mass << "\n";
                cout << "  beta = " << beta << "\n";
                cout << "  gamma = " << gamma << "\n";
                cout << "  charge = " << charge << "\n";
                cout << "  mnoc = " << mnoc << "\n";
            }
        }
        // AZ
 
        // Distance from entry point (x,y,z) to exit point (along z-axis)
        Double_t z_distance_from_in_to_out;
 
        // Condition: track direction along z-axis
        if (pz > 0.0) {
            z_distance_from_in_to_out = ModuleThickness - (z - ZStartModulePosition);
        } else {
            z_distance_from_in_to_out = z - ZStartModulePosition;
        }
 
        // Condition: track distance along z-axis must be not zero
        if (z_distance_from_in_to_out <= 0.0) {
            LOG(warn) << "BmnGemStripModule::AddRealPointFull(): Point (" << x << " : " << y << " : " << z
                      << ") has a track with an incorrect length";
            return false;
        }
 
        // Scale coefficient of the track vector (the path from entry point to exit point)
        Double_t vector_coeff = fabs(z_distance_from_in_to_out / pz);
 
        // Components of the track vector
        Double_t x_vec_from_in_to_out = vector_coeff * px;
        Double_t y_vec_from_in_to_out = vector_coeff * py;
        Double_t z_vec_from_in_to_out = vector_coeff * pz;
 
        // Exit point coordinates (x_out, y_out, z_out)
        Double_t x_out = x + x_vec_from_in_to_out;
        Double_t y_out = y + y_vec_from_in_to_out;
        // Double_t z_out = z + z_vec_from_in_to_out;
 
        // Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
        // Condition: x-coordinate of the exit point is inside the module
        if (x_out < XMinModule || x_out > XMaxModule) {
            if (x_out < XMinModule)
                x_out = XMinModule;
            if (x_out > XMaxModule)
                x_out = XMaxModule;
 
            x_vec_from_in_to_out = x_out - x;
            vector_coeff = x_vec_from_in_to_out / px;
 
            y_vec_from_in_to_out = vector_coeff * py;
            z_vec_from_in_to_out = vector_coeff * pz;
 
            y_out = y + y_vec_from_in_to_out;
            // z_out = z + z_vec_from_in_to_out;
        }
        // Condition: y-coordinate of the exit point is inside the module
        if (y_out < YMinModule || y_out > YMaxModule) {
            if (y_out < YMinModule)
                y_out = YMinModule;
            if (y_out > YMaxModule)
                y_out = YMaxModule;
 
            y_vec_from_in_to_out = y_out - y;
            vector_coeff = y_vec_from_in_to_out / py;
 
            x_vec_from_in_to_out = vector_coeff * px;
            z_vec_from_in_to_out = vector_coeff * pz;
 
            x_out = x + x_vec_from_in_to_out;
            // z_out = z + z_vec_from_in_to_out;
        }
        //----------------------------------------------------------------------
 
        // Common track length (Fix: a square root!)
        Double_t track_length =
            std::sqrt((x_vec_from_in_to_out) * (x_vec_from_in_to_out) + (y_vec_from_in_to_out) * (y_vec_from_in_to_out)
                      + (z_vec_from_in_to_out) * (z_vec_from_in_to_out));
 
        Double_t current_length = 0.0;         // traversed length (counter)
        Double_t current_length_ratio = 0.0;   // ratio of the traversed length to common track length (counter)
 
        // Int_t collisions_cnt = 0; //total collision counter
        Double_t current_step = 0.0;   // current distance between two collision points
 
        // Collection of collision points
        std::vector<CollPoint> collision_points;
 
        while (current_length < track_length) {
 
            // Reading the MCD value (mean collision distance = mean free flight path) [cm]
            MCD = BmnGemStripMedium::GetInstance().MCD;
            MCD = 1 / mnoc;   // AZ-110124
 
            current_step = gRandom->Exp(MCD);
            current_length += current_step;
 
            if (current_length > track_length)
                break;
 
            current_length_ratio = current_length / track_length;
 
            Double_t current_x = x + current_length_ratio * x_vec_from_in_to_out;
            Double_t current_y = y + current_length_ratio * y_vec_from_in_to_out;
            Double_t current_z = z + current_length_ratio * z_vec_from_in_to_out;
 
            collision_points.push_back(CollPoint(current_x, current_y, current_z));
 
            // collisions_cnt++;
        }
 
        // Each level - distance to the readout plane
        Double_t level1 = InductionGapThickness;
        Double_t level2 = InductionGapThickness + SecondTransferGapThickness;
        Double_t level3 = InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness;
        // Double_t level4 =
        // InductionGapThickness+SecondTransferGapThickness+FirstTransferGapThickness+DriftGapThickness; // not used yet
 
        // Mean electron shift along x-axis (under the influence of the Lorentz force)
        Double_t xmean;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
 
        // Sigma electron smearing
        Double_t sigma;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3+ p4*x^4 + p5*x^5)
 
        // cout << "$$$ BmnGemStripMedium::GetInstance().NameOfCurrentMedium = " <<
        // BmnGemStripMedium::GetInstance().NameOfCurrentMedium << "\n"; //test cout
 
        // Reading the medium configuration -------------------------------------
        Double_t p0_xmean = BmnGemStripMedium::GetInstance().p0_xmean;
        Double_t p1_xmean = BmnGemStripMedium::GetInstance().p1_xmean;
        Double_t p2_xmean = BmnGemStripMedium::GetInstance().p2_xmean;
        Double_t p3_xmean = BmnGemStripMedium::GetInstance().p3_xmean;
 
        Double_t p0_sigma = BmnGemStripMedium::GetInstance().p0_sigma;
        Double_t p1_sigma = BmnGemStripMedium::GetInstance().p1_sigma;
        Double_t p2_sigma = BmnGemStripMedium::GetInstance().p2_sigma;
        Double_t p3_sigma = BmnGemStripMedium::GetInstance().p3_sigma;
        Double_t p4_sigma = BmnGemStripMedium::GetInstance().p4_sigma;
        Double_t p5_sigma = BmnGemStripMedium::GetInstance().p5_sigma;
        //----------------------------------------------------------------------
 
        // Strip cluster in each strip layer
        Int_t n_strip_layers = StripLayers.size();
        vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
        // Double_t signal_coeff = signal*1E-5;
 
        // Electron avalanche producing for each collision point at the track
        for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
 
            Double_t xcoll = collision_points[icoll].x;
            Double_t ycoll = collision_points[icoll].y;
            Double_t zcoll = collision_points[icoll].z;
 
            Double_t zdist;   // current z-distance to the readout
 
            // Find z-distance to the readout depending on the electron drift direction
            if (ElectronDriftDirection == ForwardZAxisEDrift) {
                zdist = (ModuleThickness + ZStartModulePosition) - zcoll;
            } else {
                zdist = zcoll - ZStartModulePosition;
            }
            //------------------------------------------------------------------
 
            // polynomials
            xmean = p0_xmean + p1_xmean * zdist + p2_xmean * zdist * zdist + p3_xmean * zdist * zdist * zdist;
            sigma = p0_sigma + p1_sigma * zdist + p2_sigma * zdist * zdist + p3_sigma * zdist * zdist * zdist
                    + p4_sigma * zdist * zdist * zdist * zdist + p5_sigma * zdist * zdist * zdist * zdist * zdist;
 
            // Number of electrons in the current collision
            Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
            if (n_electrons_cluster < 1)
                n_electrons_cluster = 1;   // min
            if (n_electrons_cluster > 6)
                n_electrons_cluster = 6;   // max
 
            for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
 
                // Electron gain in each GEM cascade
                // Polya distribution is better, but Exponential is good too in our case
                Double_t gain_gem1 =
                    gRandom->Exp(15);   //...->Exp(V), where V is the mean value of the exponential distribution
                Double_t gain_gem2 = gRandom->Exp(15);
                Double_t gain_gem3 = gRandom->Exp(15);
 
                if (gain_gem1 < 1.0)
                    gain_gem1 = 1.0;
                if (gain_gem2 < 1.0)
                    gain_gem2 = 1.0;
                if (gain_gem3 < 1.0)
                    gain_gem3 = 1.0;
 
                int total_gain = 0;
 
                if (zdist < level1) {
                    total_gain = 1.0;
                }
                if (zdist >= level1 && zdist < level2) {
                    total_gain = gain_gem3;
                }
                if (zdist >= level2 && zdist < level3) {
                    total_gain = gain_gem3 * gain_gem2;
                }
                if (zdist >= level3) {
                    total_gain = gain_gem3 * gain_gem2 * gain_gem1;
                }
 
                // Projection of the current electron on the readout (x,y-coordinates)
                double x_readout, y_readout;
 
                for (int igain = 0; igain < total_gain; ++igain) {
 
                    // x-shift of the electon depends on the electron drift direction
                    if (ElectronDriftDirection == ForwardZAxisEDrift) {
                        x_readout = gRandom->Gaus(xcoll + xmean, sigma);
                    } else {
                        x_readout = gRandom->Gaus(xcoll - xmean, sigma);
                    }
 
                    y_readout = gRandom->Gaus(ycoll, sigma);
 
                    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                        // Condition: end electron position must be inside the layer
                        if (!StripLayers[ilayer].IsPointInsideStripLayer(x_readout, y_readout))
                            continue;
 
                        Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_readout, y_readout);
                        // AZ-100124 - account for strip width
                        Double_t dstrip = strip_pos - (int)strip_pos;
                        if (TMath::Abs(StripLayers[ilayer].GetAngleDeg()) < 1.0) {
                            // Vertical strips
                            // if (strip_pos - (int)strip_pos > 0.2) continue; //  strip width 0.16 mm
                            if (dstrip >= 0.4 && dstrip < 0.6)
                                continue;   //
                        } else {
                            // Rotated strips
                            // if (strip_pos - (int)strip_pos > 0.85) continue; //  strip width 0.68 mm
                            if (dstrip >= 0.3 && dstrip < 0.7)
                                continue;   //
                        }
                        // AZ
 
                        // Add strips with signals to lower and upper clusters
                        // StripLayers[ilayer].AddStripSignal((Int_t)strip_pos, 1.0); //Instead of 1.0 we can use Gain
                        // in future
                        cluster_layers[ilayer].AddStrip((Int_t)strip_pos, 1.0);
                        // cluster_layers[ilayer].AddStrip((Int_t)strip_pos, Gain*signal_coeff);
                    }
                    /*
                    //Electron points on the readout (testing) -------------------------------------
                                        XElectronPos.push_back(x_readout);
                                        YElectronPos.push_back(y_readout);
                                        ElectronSignal.push_back(1.0);
                    //------------------------------------------------------------------------------
                    */
                }
            }
        }
 
        // Calculate cluster parameters -----------------------------------------
        vector<Double_t> cluster_mean_position(n_strip_layers, 0.0);
        vector<Double_t> cluster_total_signal(n_strip_layers, 0.0);
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
            if (cluster_layers[ilayer].GetClusterSize() > 0) {
                for (int i = 0; i < cluster_layers[ilayer].GetClusterSize(); ++i) {
                    cluster_mean_position[ilayer] += (cluster_layers[ilayer].Strips[i] + 0.5)
                                                     * cluster_layers[ilayer].Signals[i];   // as sum of all positions
                    cluster_total_signal[ilayer] += cluster_layers[ilayer].Signals[i];
                }
                cluster_mean_position[ilayer] /= cluster_total_signal[ilayer];
            }
 
            cluster_layers[ilayer].MeanPosition = cluster_mean_position[ilayer];   // mean cluster position
            cluster_layers[ilayer].TotalSignal = cluster_total_signal[ilayer];     // total signal of the cluster
        }
 
        if (ElectronDriftDirection == ForwardZAxisEDrift) {
            if (pz > 0.0) {
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
                }
            } else {
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    cluster_layers[ilayer].OriginPosition =
                        StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
                }
            }
        } else {
            if (pz > 0.0) {
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    cluster_layers[ilayer].OriginPosition =
                        StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
                }
            } else {
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
                }
            }
        }
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            cluster_layers[ilayer].PositionResidual =
                cluster_layers[ilayer].MeanPosition - cluster_layers[ilayer].OriginPosition;
        }
        //----------------------------------------------------------------------
        /*
        //Testing ----------------------------------------------------------------------
                if(AddedClusters.size() == 0) {
                    for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                        AddedClusters.push_back(vector<StripCluster>());
                    }
                }
 
                for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
                }
        //------------------------------------------------------------------------------
        */
        // Add the correct clusters to the strip layers -------------------------
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
                Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
                Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
                StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
            }
        }
        //----------------------------------------------------------------------
 
        // Fill strip matches ---------------------------------------------------
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
                Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
                Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
                StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                        refID);
            }
        }
        //----------------------------------------------------------------------
 
        RealPointsX.push_back(x);
        RealPointsY.push_back(y);
        RealPointsXout.push_back(x_out);
        RealPointsYout.push_back(y_out);
        RealPointsMC.push_back(refID);
 
        return true;
    } else {
        LOG(warn) << "BmnGemStripModule::AddRealPointFull(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
}
// -----------------------------------------------------------------------------
 
// ANN-based cluster generation: the model returns a 2D image of the signal around the MC point
Bool_t BmnGemStripModule::AddRealPointFull_ANNGenerator(Double_t x,
                                                        Double_t y,
                                                        Double_t z,
                                                        Double_t px,
                                                        Double_t py,
                                                        Double_t pz,
                                                        Double_t signal,
                                                        Int_t refID,
                                                        Int_t pdg_code,
                                                        Double_t bfield,
                                                        Bool_t onlyStandard,
                                                        vector<Int_t> segs)
{
    Double_t maxThetaAngle = 50.0;   // max absolute angle (deg) of track inclination from beam axis
 
    // Model format: [px, py, pz, drift, bfield, s0..s15] -> image nx_bins x ny_bins in [-1, 1]
    const BmnGemStripFdeepModel* ann_model = nullptr;   // 2D fdeep model (taken below, only if the ANN is used)
    const Int_t n_track_segments = 64;                  // track segmentation (SegmentedChargeDistribution)
    const Int_t n_segs_features = 16;                   // NOTE: the first 16 of 64 segments (as before)
    const Int_t n_ann_inputs = 5 + n_segs_features;     // total number of input features
    const Int_t nx_bins = 80, ny_bins = 80;             // model output: image size
    const Double_t dx_area = 1.25, dy_area = 1.25;      // half-size of the image area (cm)
 
    ResetElectronPointsAndClusters();   // test
 
    SegmentedChargeDistribution.assign(n_track_segments, 0);
 
    // gRandom is not reseeded (as in AddRealPointFull): the seed of the run is set by the caller
 
    // Condition: a start point is inside the module and outside its dead zone (in strip layers)
    if (!IsPointInsideModule(x, y)) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
 
    TVector3 mom(px, py, pz);
    TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg_code);
 
    Double_t mnoc = 1.0, mass = 3.7283999;   // He4
    if (part)
        mass = part->Mass();
    if (mass > 0.0001) {
        TLorentzVector lorv;
        lorv.SetVectM(mom, mass);
        Double_t beta = lorv.Beta();
        Double_t gamma = lorv.Gamma();
        Double_t charge = (part) ? part->Charge() / 3 : 1;
        mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
        if (mnoc < 0.5)
            return kFALSE;
        if (MaxClustersPerCm > 0.0 && mnoc > MaxClustersPerCm)
            mnoc = MaxClustersPerCm;
    }
 
    // Distance from entry point (x,y,z) to exit point (along z-axis)
    Double_t z_distance_from_in_to_out;
 
    // Condition: track direction along z-axis
    if (pz > 0.0) {
        z_distance_from_in_to_out = ModuleThickness - (z - ZStartModulePosition);
    } else {
        z_distance_from_in_to_out = z - ZStartModulePosition;
    }
 
    // Condition: track distance along z-axis must be not zero
    if (z_distance_from_in_to_out <= 0.0) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") has a track with an incorrect length";
        return false;
    }
 
    // Scale coefficient of the track vector (the path from entry point to exit point)
    Double_t vector_coeff = fabs(z_distance_from_in_to_out / pz);
 
    // Components of the track vector
    Double_t x_vec_from_in_to_out = vector_coeff * px;
    Double_t y_vec_from_in_to_out = vector_coeff * py;
    Double_t z_vec_from_in_to_out = vector_coeff * pz;
 
    // Exit point coordinates (x_out, y_out, z_out)
    Double_t x_out = x + x_vec_from_in_to_out;
    Double_t y_out = y + y_vec_from_in_to_out;
    Double_t z_out = z + z_vec_from_in_to_out;
 
    // Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
    // Condition: x-coordinate of the exit point is inside the module
    if (x_out < XMinModule || x_out > XMaxModule) {
        if (x_out < XMinModule)
            x_out = XMinModule;
        if (x_out > XMaxModule)
            x_out = XMaxModule;
 
        x_vec_from_in_to_out = x_out - x;
        vector_coeff = x_vec_from_in_to_out / px;
 
        y_vec_from_in_to_out = vector_coeff * py;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        y_out = y + y_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    // Condition: y-coordinate of the exit point is inside the module
    if (y_out < YMinModule || y_out > YMaxModule) {
        if (y_out < YMinModule)
            y_out = YMinModule;
        if (y_out > YMaxModule)
            y_out = YMaxModule;
 
        y_vec_from_in_to_out = y_out - y;
        vector_coeff = y_vec_from_in_to_out / py;
 
        x_vec_from_in_to_out = vector_coeff * px;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        x_out = x + x_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    //----------------------------------------------------------------------
 
    // Common track length
    Double_t track_length =
        std::sqrt((x_vec_from_in_to_out) * (x_vec_from_in_to_out) + (y_vec_from_in_to_out) * (y_vec_from_in_to_out)
                  + (z_vec_from_in_to_out) * (z_vec_from_in_to_out));
 
    Double_t current_length = 0.0;         // traversed length (counter)
    Double_t current_length_ratio = 0.0;   // ratio of the traversed length to common track length (counter)
 
    Int_t collisions_cnt = 0;      // total collision counter
    Double_t current_step = 0.0;   // current distance between two collision points
 
    // Collection of collision points
    std::vector<CollPoint> collision_points;
 
    while (current_length < track_length) {
 
        // Reading the MCD value (mean collision distance = mean free flight path) [cm]
        MCD = BmnGemStripMedium::GetInstance().MCD;
        MCD = 1 / mnoc;   // AZ-110124
 
        current_step = gRandom->Exp(MCD);
        current_length += current_step;
 
        if (current_length > track_length)
            break;
 
        current_length_ratio = current_length / track_length;
 
        Double_t current_x = x + current_length_ratio * x_vec_from_in_to_out;
        Double_t current_y = y + current_length_ratio * y_vec_from_in_to_out;
        Double_t current_z = z + current_length_ratio * z_vec_from_in_to_out;
 
        collision_points.push_back(CollPoint(current_x, current_y, current_z));
 
        collisions_cnt++;
    }
 
    // Each level - distance to the readout plane
    Double_t level1 = InductionGapThickness;
    Double_t level2 = InductionGapThickness + SecondTransferGapThickness;
    Double_t level3 = InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness;
 
    // Mean electron shift along x-axis (under the influence of the Lorentz force)
    Double_t xmean;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
 
    // Sigma electron smearing
    Double_t sigma;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3+ p4*x^4 + p5*x^5)
 
    // Reading the medium configuration -------------------------------------
    Double_t p0_xmean = BmnGemStripMedium::GetInstance().p0_xmean;
    Double_t p1_xmean = BmnGemStripMedium::GetInstance().p1_xmean;
    Double_t p2_xmean = BmnGemStripMedium::GetInstance().p2_xmean;
    Double_t p3_xmean = BmnGemStripMedium::GetInstance().p3_xmean;
 
    Double_t p0_sigma = BmnGemStripMedium::GetInstance().p0_sigma;
    Double_t p1_sigma = BmnGemStripMedium::GetInstance().p1_sigma;
    Double_t p2_sigma = BmnGemStripMedium::GetInstance().p2_sigma;
    Double_t p3_sigma = BmnGemStripMedium::GetInstance().p3_sigma;
    Double_t p4_sigma = BmnGemStripMedium::GetInstance().p4_sigma;
    Double_t p5_sigma = BmnGemStripMedium::GetInstance().p5_sigma;
    //----------------------------------------------------------------------
 
    Double_t r_vec = TMath::Sqrt(px * px + py * py + pz * pz);
    Double_t theta_deg_abs = TMath::ACos(TMath::Abs(pz) / r_vec) * TMath::RadToDeg();
 
    // fdeep model: loaded on the first request (not needed with onlyStandard), shared between modules;
    // its input / output sizes are checked once for each model
    std::shared_ptr<const BmnGemStripFdeepModel> ann_holder;
    if (!onlyStandard) {
        ann_holder = BmnGemStripFdeepModel::GetShared();
        ann_model = ann_holder.get();
 
        static std::mutex checked_mutex;
        static std::map<const BmnGemStripFdeepModel*,
                        std::pair<std::shared_ptr<const BmnGemStripFdeepModel>, std::string>>
            checked;   // model -> (holder, why it is not suitable: empty if suitable)
 
        std::lock_guard<std::mutex> lock(checked_mutex);
        if (!ann_model) {
            LOG(warn) << "BmnGemStripModule::" << __func__
                      << "(): fdeep model is not available, the standard method is used";
        } else {
            auto it = checked.find(ann_model);
            if (it == checked.end()) {
                const Int_t expected_out = nx_bins * ny_bins;
                std::string why;
                if ((ann_model->GetInputDim() != n_ann_inputs) || (ann_model->GetOutputDim() != expected_out))
                    why = "input " + std::to_string(ann_model->GetInputDim()) + ", expected "
                          + std::to_string(n_ann_inputs) + "; output " + std::to_string(ann_model->GetOutputDim())
                          + ", expected " + std::to_string(expected_out);
                else
                    LOG(info) << "BmnGemStripModule::" << __func__ << "(): fdeep model " << ann_model->GetFilePath()
                              << " is used";
                it = checked.emplace(ann_model, std::make_pair(ann_holder, why)).first;
            }
            if (!it->second.second.empty()) {
                LOG(warn) << "BmnGemStripModule::" << __func__ << "(): the model " << ann_model->GetFilePath()
                          << " is not suitable (" << it->second.second << "), the standard method is used";
                ann_model = nullptr;
            }
        }
    }
 
    // Which method is used for this point: standard (electrons -> strips) or ANN (model -> strips)
    // (without a loaded model - e.g. built without fdeep - the standard method is used)
    const Bool_t useANN = (theta_deg_abs < maxThetaAngle) && !onlyStandard && ann_model;
 
    // Initialization of a strip cluster in each strip layer
    Int_t n_strip_layers = StripLayers.size();   // number of layers in the module
    vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
    // Electron avalanche producing for each collision point at the track
    // (in the ANN mode only the segmented charge distribution is filled here)
    for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
 
        Double_t xcoll = collision_points[icoll].x;
        Double_t ycoll = collision_points[icoll].y;
        Double_t zcoll = collision_points[icoll].z;
 
        Double_t zdist;   // current z-distance to the readout
 
        // Find z-distance to the readout depending on the electron drift direction
        if (ElectronDriftDirection == ForwardZAxisEDrift) {
            zdist = (ModuleThickness + ZStartModulePosition) - zcoll;
        } else {
            zdist = zcoll - ZStartModulePosition;
        }
        //------------------------------------------------------------------
 
        // polynomials
        xmean = p0_xmean + p1_xmean * zdist + p2_xmean * zdist * zdist + p3_xmean * zdist * zdist * zdist;
        sigma = p0_sigma + p1_sigma * zdist + p2_sigma * zdist * zdist + p3_sigma * zdist * zdist * zdist
                + p4_sigma * zdist * zdist * zdist * zdist + p5_sigma * zdist * zdist * zdist * zdist * zdist;
 
        // Number of electrons in the current collision
        Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
        if (n_electrons_cluster < 1)
            n_electrons_cluster = 1;   // min
        if (n_electrons_cluster > 6)
            n_electrons_cluster = 6;   // max
 
        // Filling charge distribution -------------------------------------
        Double_t dist_from_start = 0.0;
        if (pz > 0.0) {
            dist_from_start =
                std::sqrt((xcoll - x) * (xcoll - x) + (ycoll - y) * (ycoll - y) + (zcoll - z) * (zcoll - z));
        } else {
            dist_from_start = std::sqrt((xcoll - x_out) * (xcoll - x_out) + (ycoll - y_out) * (ycoll - y_out)
                                        + (zcoll - z_out) * (zcoll - z_out));
        }
        Int_t segment = (Int_t)(dist_from_start / track_length * n_track_segments);
        if (segment >= n_track_segments)
            segment = n_track_segments - 1;
        //------------------------------------------------------------------
 
        for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
 
            // Electron gain in each GEM cascade
            Double_t gain_gem1 = gRandom->Exp(15);
            Double_t gain_gem2 = gRandom->Exp(15);
            Double_t gain_gem3 = gRandom->Exp(15);
 
            if (gain_gem1 < 1.0)
                gain_gem1 = 1.0;
            if (gain_gem2 < 1.0)
                gain_gem2 = 1.0;
            if (gain_gem3 < 1.0)
                gain_gem3 = 1.0;
 
            int total_gain = 0;
 
            if (zdist < level1) {
                total_gain = 1.0;
            }
            if (zdist >= level1 && zdist < level2) {
                total_gain = gain_gem3;
            }
            if (zdist >= level2 && zdist < level3) {
                total_gain = gain_gem3 * gain_gem2;
            }
            if (zdist >= level3) {
                total_gain = gain_gem3 * gain_gem2 * gain_gem1;
            }
 
            SegmentedChargeDistribution[segment] += total_gain;
 
            // ANN mode: electrons are not projected on the readout, strips are taken from the model
            if (useANN)
                continue;
 
            // Standard method ---------------------------------------------------
            // Projection of the current electron on the readout (x,y-coordinates)
            double x_readout, y_readout;
 
            for (int igain = 0; igain < total_gain; ++igain) {
 
                // x-shift of the electon depends on the electron drift direction
                if (ElectronDriftDirection == ForwardZAxisEDrift) {
                    x_readout = gRandom->Gaus(xcoll + xmean, sigma);
                } else {
                    x_readout = gRandom->Gaus(xcoll - xmean, sigma);
                }
 
                y_readout = gRandom->Gaus(ycoll, sigma);
 
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    // Condition: end electron position must be inside the layer
                    if (!StripLayers[ilayer].IsPointInsideStripLayer(x_readout, y_readout))
                        continue;
 
                    Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_readout, y_readout);
                    // AZ-100124 - account for strip width
                    Double_t dstrip = strip_pos - (int)strip_pos;
                    if (TMath::Abs(StripLayers[ilayer].GetAngleDeg()) < 1.0) {
                        // Vertical strips
                        if (dstrip >= 0.4 && dstrip < 0.6)
                            continue;
                    } else {
                        // Rotated strips
                        if (dstrip >= 0.3 && dstrip < 0.7)
                            continue;
                    }
                    // AZ
 
                    cluster_layers[ilayer].AddStrip((Int_t)strip_pos, 1.0);
                }
 
                // Electron points on the readout (testing) -------------
                XElectronPos.push_back(x_readout);
                YElectronPos.push_back(y_readout);
                ElectronSignal.push_back(1.0);
                //------------------------------------------------------
            }
        }
    }
    //----------------------------------------------------------------------
 
    // ANN generation ------------------------------------------------------
    if (useANN) {
 
        // Normalization (px,py,pz), where r_vec must be equal to 1
        Double_t px_norm = px / r_vec;
        Double_t py_norm = py / r_vec;
        Double_t pz_norm = pz / r_vec;
 
        // Inversion for z-reverse direction particles: the same track segment passed from (x_out, y_out)
        if (pz < 0.0) {
            px_norm = -px_norm;
            py_norm = -py_norm;
            pz_norm = -pz_norm;
        }
 
        Double_t drift_direction = 0;   // 0 is ForwardZAxisEDrift
        if (ElectronDriftDirection == BackwardZAxisEDrift)
            drift_direction = 1;
 
        // Set seg. vector explicitly from parameter segs if it is not empty
        if (!segs.empty()) {
            SegmentedChargeDistribution = segs;
        }
 
        if (SegmentedChargeDistribution.size() < static_cast<size_t>(n_segs_features)) {
            LOG(error) << "BmnGemStripModule::" << __func__
                       << "(): segment vector size = " << SegmentedChargeDistribution.size() << ", the model requires "
                       << n_segs_features;
            return false;
        }
 
        // Normalize seg. vector: [0, max] -> [0, 1] (max over the whole vector)
        Double_t max_segment_val =
            *(std::max_element(SegmentedChargeDistribution.begin(), SegmentedChargeDistribution.end()));
        vector<Float_t> segs_norm(n_segs_features, 0.0f);
        if (max_segment_val > 0.0) {
            for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
                segs_norm[iseg] = static_cast<Float_t>(SegmentedChargeDistribution[iseg] / max_segment_val);
        }
 
        // Total charge in segment vector
        Double_t total_charge = 0.0;
        for (size_t i = 0; i < SegmentedChargeDistribution.size(); ++i)
            total_charge += SegmentedChargeDistribution[i];
 
        // MC point: the track start point along the (inverted) momentum
        Double_t x_mc = x;
        Double_t y_mc = y;
        if (pz < 0.0) {
            x_mc = x_out;
            y_mc = y_out;
        }
 
        // Common part of the model input: [px, py, pz, drift, bfield, s0..s{n_segs_features-1}]
        std::vector<float> input_data;
        input_data.reserve(n_ann_inputs);
        input_data.push_back(static_cast<float>(px_norm));
        input_data.push_back(static_cast<float>(py_norm));
        input_data.push_back(static_cast<float>(pz_norm));
        input_data.push_back(static_cast<float>(drift_direction));
        input_data.push_back(static_cast<float>(bfield));
        input_data.insert(input_data.end(), segs_norm.begin(), segs_norm.end());
 
        if (total_charge > 0.0) {
            // 2D model: image of the signal around the MC point -> random points in its bins -> strips
            const auto img = ann_model->Predict(input_data);
            if (img.size() != static_cast<std::size_t>(nx_bins * ny_bins)) {
                LOG(error) << "BmnGemStripModule::" << __func__ << "(): model output size = " << img.size()
                           << ", expected " << nx_bins * ny_bins;
            } else {
                Double_t x_bin_size = 2 * dx_area / nx_bins;
                Double_t y_bin_size = 2 * dy_area / ny_bins;
                Double_t x_min_area = x_mc - dx_area;
                Double_t y_min_area = y_mc - dy_area;
 
                // Normalize content of img: [-1,1] -> [0,1]
                std::vector<Double_t> norm_content(nx_bins * ny_bins, 0.0);
                Double_t sum_norm = 0.0;
                for (Int_t idx = 0; idx < nx_bins * ny_bins; ++idx) {
                    norm_content[idx] = (img.at(idx) - (-1.0)) / (1.0 - (-1.0));
                    sum_norm += norm_content[idx];
                }
 
                // Scale coeff for total signal distribution
                const Double_t ann_correction = 0.8;   // Set to 1.0 to disable correction
                Double_t signal_scale = 0.0;
                if (sum_norm > 0.0)
                    signal_scale = (total_charge / sum_norm) * ann_correction;
 
                Double_t threshold = 0.0001 * total_charge;
 
                // Image rows go from the top (iy = ny_bins - 1) to the bottom
                Int_t idx = 0;
                for (Int_t iy = ny_bins - 1; iy >= 0; --iy) {
                    for (Int_t ix = 0; ix < nx_bins; ++ix, ++idx) {
                        Double_t content = norm_content[idx] * signal_scale;
                        if (content <= threshold)
                            continue;
 
                        Double_t x_bin_coord = ix * x_bin_size + x_bin_size * 0.5 + x_min_area;
                        Double_t y_bin_coord = iy * y_bin_size + y_bin_size * 0.5 + y_min_area;
 
                        Int_t n_rand_points = 10;
                        Double_t distr_content = content / n_rand_points;
                        for (Int_t irp = 0; irp < n_rand_points; ++irp) {
                            Double_t x_distr =
                                gRandom->Uniform(x_bin_coord - x_bin_size * 0.5, x_bin_coord + x_bin_size * 0.5);
                            Double_t y_distr =
                                gRandom->Uniform(y_bin_coord - y_bin_size * 0.5, y_bin_coord + y_bin_size * 0.5);
 
                            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                                if (!StripLayers[ilayer].IsPointInsideStripLayer(x_distr, y_distr))
                                    continue;
                                Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_distr, y_distr);
                                cluster_layers[ilayer].AddStrip((Int_t)strip_pos, distr_content);
                            }
 
                            XElectronPos.push_back(x_distr);
                            YElectronPos.push_back(y_distr);
                            ElectronSignal.push_back(distr_content);
                        }
                    }
                }
            }
        }
    }
    //----------------------------------------------------------------------
 
    // Calculate cluster parameters --------------------------------------------
    vector<Double_t> cluster_mean_position(n_strip_layers, 0.0);
    vector<Double_t> cluster_total_signal(n_strip_layers, 0.0);
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
        if (cluster_layers[ilayer].GetClusterSize() > 0) {
 
            for (int i = 0; i < cluster_layers[ilayer].GetClusterSize(); ++i) {
                cluster_mean_position[ilayer] += (cluster_layers[ilayer].Strips[i] + 0.5)
                                                 * cluster_layers[ilayer].Signals[i];   // as sum of all positions
                cluster_total_signal[ilayer] += cluster_layers[ilayer].Signals[i];
            }
            cluster_mean_position[ilayer] /= cluster_total_signal[ilayer];
        }
 
        cluster_layers[ilayer].MeanPosition = cluster_mean_position[ilayer];   // mean cluster position
        cluster_layers[ilayer].TotalSignal = cluster_total_signal[ilayer];     // total signal of the cluster
    }
 
    if (ElectronDriftDirection == ForwardZAxisEDrift) {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        }
    } else {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        cluster_layers[ilayer].PositionResidual =
            cluster_layers[ilayer].MeanPosition - cluster_layers[ilayer].OriginPosition;
        cluster_layers[ilayer].SetType(StripLayers[ilayer].GetType());
    }
    //----------------------------------------------------------------------
 
    // Testing ----------------------------------------------------------------------
    if (AddedClusters.size() == 0) {
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            AddedClusters.push_back(vector<StripCluster>());
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
    }
    //------------------------------------------------------------------------------
 
    // Add the correct clusters to the strip layers -------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
        }
    }
    //----------------------------------------------------------------------
 
    // Fill strip matches ---------------------------------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                    refID);
        }
    }
    //----------------------------------------------------------------------
 
    RealPointsX.push_back(x);
    RealPointsY.push_back(y);
    RealPointsXout.push_back(x_out);
    RealPointsYout.push_back(y_out);
    RealPointsMC.push_back(refID);
 
    return true;
}
//------------------------------------------------------------------------------
 
// Strip generation with the fdeep model (BmnGemStripFdeepModel): standard method or ANN (|theta| < kANNMaxThetaAngle).
// Generation data only with SetFillGenerationData(kTRUE);
Bool_t BmnGemStripModule::AddRealPointFull_ANNStripGenerator_Fdeep(Double_t x,
                                                                   Double_t y,
                                                                   Double_t z,
                                                                   Double_t px,
                                                                   Double_t py,
                                                                   Double_t pz,
                                                                   Double_t signal,
                                                                   Int_t refID,
                                                                   Int_t pdg_code,
                                                                   Double_t bfield,
                                                                   Bool_t onlyStandard,
                                                                   vector<Int_t> segs)
{
    const Double_t maxThetaAngle = kANNMaxThetaAngle;   // max |theta| (deg) of the track for the ANN
    const Bool_t kParallelANNLayers = false;   // [PARALLEL-ANN] true - ANN prediction for the strip layers in parallel
 
    // Model format: [px, py, pz, drift, bfield, s0..s63, strip_layer, strip_order, center_frac] -> window of N
    // strips (N odd, from the model) centered at the middle of the drift gap (GetANNWindowCenterPoint)
    const Int_t n_track_segments = 64;                    // track segmentation (segmented charge distribution)
    const Int_t n_segs_features = n_track_segments;       // segments passed to the model
    const Int_t n_ann_inputs = 5 + n_segs_features + 3;   // total number of input features
    Int_t n_ann_strips = 0;   // model output: strip window size (odd), taken from the model below
 
    // Test outputs are cleared always (cheap if they are empty)
    ResetElectronPointsAndClusters();
 
    // fdeep model: shared between modules, checked once for each model -------
    std::shared_ptr<const BmnGemStripFdeepModel> ann_holder;
    const BmnGemStripFdeepModel* ann_model = nullptr;
    if (!onlyStandard) {
        ann_holder = BmnGemStripFdeepModel::GetShared();
        ann_model = ann_holder.get();
 
        static std::mutex checked_mutex;
        static std::map<const BmnGemStripFdeepModel*,
                        std::pair<std::shared_ptr<const BmnGemStripFdeepModel>, std::string>>
            checked;   // model -> (holder, why it is not suitable: empty if suitable)
 
        std::lock_guard<std::mutex> lock(checked_mutex);
        if (!ann_model) {
            LOG(warn) << "BmnGemStripModule::" << __func__
                      << "(): fdeep model is not available, the standard method is used";
        } else {
            auto it = checked.find(ann_model);
            if (it == checked.end()) {
                std::string why;
                if (ann_model->GetInputDim() != n_ann_inputs)
                    why = "input " + std::to_string(ann_model->GetInputDim()) + ", expected "
                          + std::to_string(n_ann_inputs);
                else
                    why = CheckANNStripWindowSize(ann_model->GetOutputDim());
                if (why.empty())
                    LOG(info) << "BmnGemStripModule::" << __func__ << "(): fdeep model " << ann_model->GetFilePath()
                              << " is used (ANN for |theta| < " << maxThetaAngle << " deg, window "
                              << ann_model->GetOutputDim() << " strips)";
                it = checked.emplace(ann_model, std::make_pair(ann_holder, why)).first;
            }
            if (!it->second.second.empty()) {
                LOG(warn) << "BmnGemStripModule::" << __func__ << "(): the model " << ann_model->GetFilePath()
                          << " is not suitable (" << it->second.second << "), the standard method is used";
                ann_model = nullptr;
            }
        }
    }
    if (ann_model)
        n_ann_strips = ann_model->GetOutputDim();   // window size of the model (checked above: odd, >= 3)
    //--------------------------------------------------------------------------
 
    // Condition: a start point is inside the module and outside its dead zone (in strip layers)
    if (!IsPointInsideModule(x, y)) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
 
    // Particle: pdg_code = 0 -> from the MC track of the point refID (as in AddRealPointFull) --
    Int_t pdg = pdg_code;
    if (pdg == 0) {
        FairRunSim* run = FairRunSim::Instance();
        BmnGemStripDigitizer* digi = run ? (BmnGemStripDigitizer*)run->GetTask("GemDigitizer") : nullptr;
        TClonesArray* points = digi ? digi->GetMCPoints() : nullptr;
        TClonesArray* mctracks = digi ? digi->GetMCTracks() : nullptr;
 
        CbmMCTrack* mcTr = nullptr;
        if (points && mctracks && refID >= 0 && refID < points->GetEntriesFast()) {
            CbmStsPoint* point = (CbmStsPoint*)points->UncheckedAt(refID);
            if (point && point->GetTrackID() >= 0 && point->GetTrackID() < mctracks->GetEntriesFast())
                mcTr = (CbmMCTrack*)mctracks->UncheckedAt(point->GetTrackID());
        }
 
        if (mcTr) {
            pdg = mcTr->GetPdgCode();
        } else {
            static std::atomic<bool> no_mc_reported{false};
            if (!no_mc_reported.exchange(true))
                LOG(warn) << "BmnGemStripModule::" << __func__
                          << "(): pdg_code = 0 and the MC track is not available (no GemDigitizer or MC point), "
                             "the He4 mass is used (this warning is printed once)";
        }
    }
 
    TVector3 mom(px, py, pz);
    TParticlePDG* part = (pdg != 0) ? TDatabasePDG::Instance()->GetParticle(pdg) : nullptr;
 
    Double_t mnoc = 1.0, mass = 3.7283999;   // He4
    if (part)
        mass = part->Mass();
    if (mass > 0.0001) {
        TLorentzVector lorv;
        lorv.SetVectM(mom, mass);
        Double_t beta = lorv.Beta();
        Double_t gamma = lorv.Gamma();
        Double_t charge = (part) ? part->Charge() / 3 : 1;
        mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
        if (mnoc < 0.5)
            return kFALSE;
        if (MaxClustersPerCm > 0.0 && mnoc > MaxClustersPerCm)
            mnoc = MaxClustersPerCm;
    }
    //--------------------------------------------------------------------------
 
    // Distance from entry point (x,y,z) to exit point (along z-axis)
    Double_t z_distance_from_in_to_out;
 
    // Condition: track direction along z-axis
    if (pz > 0.0) {
        z_distance_from_in_to_out = ModuleThickness - (z - ZStartModulePosition);
    } else {
        z_distance_from_in_to_out = z - ZStartModulePosition;
    }
 
    // Condition: track distance along z-axis must be not zero
    if (z_distance_from_in_to_out <= 0.0) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") has a track with an incorrect length";
        return false;
    }
 
    // Scale coefficient of the track vector (the path from entry point to exit point)
    Double_t vector_coeff = fabs(z_distance_from_in_to_out / pz);
 
    // Components of the track vector
    Double_t x_vec_from_in_to_out = vector_coeff * px;
    Double_t y_vec_from_in_to_out = vector_coeff * py;
    Double_t z_vec_from_in_to_out = vector_coeff * pz;
 
    // Exit point coordinates (x_out, y_out, z_out)
    Double_t x_out = x + x_vec_from_in_to_out;
    Double_t y_out = y + y_vec_from_in_to_out;
    Double_t z_out = z + z_vec_from_in_to_out;
 
    // Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
    // Condition: x-coordinate of the exit point is inside the module
    if (x_out < XMinModule || x_out > XMaxModule) {
        if (x_out < XMinModule)
            x_out = XMinModule;
        if (x_out > XMaxModule)
            x_out = XMaxModule;
 
        x_vec_from_in_to_out = x_out - x;
        vector_coeff = x_vec_from_in_to_out / px;
 
        y_vec_from_in_to_out = vector_coeff * py;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        y_out = y + y_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    // Condition: y-coordinate of the exit point is inside the module
    if (y_out < YMinModule || y_out > YMaxModule) {
        if (y_out < YMinModule)
            y_out = YMinModule;
        if (y_out > YMaxModule)
            y_out = YMaxModule;
 
        y_vec_from_in_to_out = y_out - y;
        vector_coeff = y_vec_from_in_to_out / py;
 
        x_vec_from_in_to_out = vector_coeff * px;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        x_out = x + x_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    //--------------------------------------------------------------------------
 
    // Center point of the ANN strip window: the track at the middle of the drift gap (computed always, also with
    // onlyStandard: it is stored with the point, see GetRealPointXcenter(), for the training samples)
    Double_t x_c = x, y_c = y;
    GetANNWindowCenterPoint(x, y, z, x_out, y_out, z_out, x_c, y_c);
 
    // Common track length
    Double_t track_length =
        std::sqrt((x_vec_from_in_to_out) * (x_vec_from_in_to_out) + (y_vec_from_in_to_out) * (y_vec_from_in_to_out)
                  + (z_vec_from_in_to_out) * (z_vec_from_in_to_out));
 
    // Which method is used for this point: standard (electrons -> strips) or ANN (model -> strips)
    // (without a suitable model or with onlyStandard the standard method is used)
    Double_t r_vec = TMath::Sqrt(px * px + py * py + pz * pz);
    Double_t theta_deg_abs = TMath::ACos(TMath::Abs(pz) / r_vec) * TMath::RadToDeg();
    const Bool_t useANN = (theta_deg_abs < maxThetaAngle) && ann_model;
 
    // ANN with external segments: the avalanche is not needed
    const Bool_t use_external_segs = useANN && !segs.empty();
    if (use_external_segs && segs.size() < static_cast<size_t>(n_segs_features)) {
        LOG(error) << "BmnGemStripModule::" << __func__ << "(): segment vector size = " << segs.size()
                   << ", the model requires " << n_segs_features;
        return false;
    }
 
    // Each level - distance to the readout plane
    Double_t level1 = InductionGapThickness;
    Double_t level2 = InductionGapThickness + SecondTransferGapThickness;
    Double_t level3 = InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness;
 
    // Mean electron shift along x-axis (under the influence of the Lorentz force)
    Double_t xmean;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
 
    // Sigma electron smearing
    Double_t sigma;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3+ p4*x^4 + p5*x^5)
 
    // Reading the medium configuration -------------------------------------
    Double_t p0_xmean = BmnGemStripMedium::GetInstance().p0_xmean;
    Double_t p1_xmean = BmnGemStripMedium::GetInstance().p1_xmean;
    Double_t p2_xmean = BmnGemStripMedium::GetInstance().p2_xmean;
    Double_t p3_xmean = BmnGemStripMedium::GetInstance().p3_xmean;
 
    Double_t p0_sigma = BmnGemStripMedium::GetInstance().p0_sigma;
    Double_t p1_sigma = BmnGemStripMedium::GetInstance().p1_sigma;
    Double_t p2_sigma = BmnGemStripMedium::GetInstance().p2_sigma;
    Double_t p3_sigma = BmnGemStripMedium::GetInstance().p3_sigma;
    Double_t p4_sigma = BmnGemStripMedium::GetInstance().p4_sigma;
    Double_t p5_sigma = BmnGemStripMedium::GetInstance().p5_sigma;
    //----------------------------------------------------------------------
 
    // Initialization of a strip cluster in each strip layer
    Int_t n_strip_layers = StripLayers.size();   // number of layers in the module
    vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
    // Segmented charge distribution of this point (always computed when the avalanche is simulated)
    vector<Int_t> seg_local(n_track_segments, 0);
 
    // Standard method: per-layer strip accumulators -------------------------
    // Signals are accumulated per strip and added to the cluster with one AddStrip() per strip:
    // AddStrip() sums the signals and keeps the strips in ascending order, the signals are integer
    // numbers of electrons -> the cluster is bitwise the same as with one AddStrip() per electron.
    vector<vector<Double_t>> strip_acc;
    vector<Int_t> acc_min, acc_max;   // range of the touched strips in each layer
    vector<char> layer_vertical;      // layer invariant: vertical strips (|angle| < 1 deg)
    if (!useANN) {
        strip_acc.resize(n_strip_layers);
        acc_min.assign(n_strip_layers, std::numeric_limits<Int_t>::max());
        acc_max.assign(n_strip_layers, -1);
        layer_vertical.assign(n_strip_layers, 0);
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            strip_acc[ilayer].assign(std::max(0, StripLayers[ilayer].GetNStrips()), 0.0);
            layer_vertical[ilayer] = TMath::Abs(StripLayers[ilayer].GetAngleDeg()) < 1.0;
        }
    }
    //--------------------------------------------------------------------------
 
    if (!use_external_segs) {
 
        Double_t current_length = 0.0;         // traversed length (counter)
        Double_t current_length_ratio = 0.0;   // ratio of the traversed length to common track length (counter)
        Double_t current_step = 0.0;           // current distance between two collision points
 
        // Collection of collision points
        std::vector<CollPoint> collision_points;
 
        // Reading the MCD value (mean collision distance = mean free flight path) [cm]
        MCD = BmnGemStripMedium::GetInstance().MCD;
        MCD = 1 / mnoc;   // AZ-110124
 
        while (current_length < track_length) {
 
            current_step = gRandom->Exp(MCD);
            current_length += current_step;
 
            if (current_length > track_length)
                break;
 
            current_length_ratio = current_length / track_length;
 
            Double_t current_x = x + current_length_ratio * x_vec_from_in_to_out;
            Double_t current_y = y + current_length_ratio * y_vec_from_in_to_out;
            Double_t current_z = z + current_length_ratio * z_vec_from_in_to_out;
 
            collision_points.push_back(CollPoint(current_x, current_y, current_z));
        }
 
        // Electron avalanche producing for each collision point at the track
        // (in the ANN mode only the segmented charge distribution is filled here)
        for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
 
            Double_t xcoll = collision_points[icoll].x;
            Double_t ycoll = collision_points[icoll].y;
            Double_t zcoll = collision_points[icoll].z;
 
            Double_t zdist;   // current z-distance to the readout
 
            // Find z-distance to the readout depending on the electron drift direction
            if (ElectronDriftDirection == ForwardZAxisEDrift) {
                zdist = (ModuleThickness + ZStartModulePosition) - zcoll;
            } else {
                zdist = zcoll - ZStartModulePosition;
            }
            //------------------------------------------------------------------
 
            // polynomials
            xmean = p0_xmean + p1_xmean * zdist + p2_xmean * zdist * zdist + p3_xmean * zdist * zdist * zdist;
            sigma = p0_sigma + p1_sigma * zdist + p2_sigma * zdist * zdist + p3_sigma * zdist * zdist * zdist
                    + p4_sigma * zdist * zdist * zdist * zdist + p5_sigma * zdist * zdist * zdist * zdist * zdist;
 
            // Number of electrons in the current collision
            Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
            if (n_electrons_cluster < 1)
                n_electrons_cluster = 1;   // min
            if (n_electrons_cluster > 6)
                n_electrons_cluster = 6;   // max
 
            // Filling charge distribution -------------------------------------
            Double_t dist_from_start = 0.0;
            if (pz > 0.0) {
                dist_from_start =
                    std::sqrt((xcoll - x) * (xcoll - x) + (ycoll - y) * (ycoll - y) + (zcoll - z) * (zcoll - z));
            } else {
                dist_from_start = std::sqrt((xcoll - x_out) * (xcoll - x_out) + (ycoll - y_out) * (ycoll - y_out)
                                            + (zcoll - z_out) * (zcoll - z_out));
            }
            Int_t segment = (Int_t)(dist_from_start / track_length * n_track_segments);
            if (segment >= n_track_segments)
                segment = n_track_segments - 1;
            //------------------------------------------------------------------
 
            for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
 
                // Electron gain in each GEM cascade (three draws always, as in the other methods)
                Double_t gain_gem1 = gRandom->Exp(15);
                Double_t gain_gem2 = gRandom->Exp(15);
                Double_t gain_gem3 = gRandom->Exp(15);
 
                if (gain_gem1 < 1.0)
                    gain_gem1 = 1.0;
                if (gain_gem2 < 1.0)
                    gain_gem2 = 1.0;
                if (gain_gem3 < 1.0)
                    gain_gem3 = 1.0;
 
                int total_gain = 0;
 
                if (zdist < level1) {
                    total_gain = 1.0;
                }
                if (zdist >= level1 && zdist < level2) {
                    total_gain = gain_gem3;
                }
                if (zdist >= level2 && zdist < level3) {
                    total_gain = gain_gem3 * gain_gem2;
                }
                if (zdist >= level3) {
                    total_gain = gain_gem3 * gain_gem2 * gain_gem1;
                }
 
                seg_local[segment] += total_gain;
 
                // ANN mode: electrons are not projected on the readout, strips are taken from the model
                if (useANN)
                    continue;
 
                // Standard method ---------------------------------------------
                // Projection of the current electron on the readout (x,y-coordinates)
                double x_readout, y_readout;
 
                for (int igain = 0; igain < total_gain; ++igain) {
 
                    // x-shift of the electon depends on the electron drift direction
                    if (ElectronDriftDirection == ForwardZAxisEDrift) {
                        x_readout = gRandom->Gaus(xcoll + xmean, sigma);
                    } else {
                        x_readout = gRandom->Gaus(xcoll - xmean, sigma);
                    }
 
                    y_readout = gRandom->Gaus(ycoll, sigma);
 
                    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                        // Condition: end electron position must be inside the layer
                        if (!StripLayers[ilayer].IsPointInsideStripLayer(x_readout, y_readout))
                            continue;
 
                        Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_readout, y_readout);
                        // AZ-100124 - account for strip width
                        Double_t dstrip = strip_pos - (int)strip_pos;
                        if (layer_vertical[ilayer]) {
                            // Vertical strips
                            if (dstrip >= 0.4 && dstrip < 0.6)
                                continue;
                        } else {
                            // Rotated strips
                            if (dstrip >= 0.3 && dstrip < 0.7)
                                continue;
                        }
                        // AZ
 
                        const Int_t strip_num = (Int_t)strip_pos;
                        if (strip_num >= 0 && strip_num < (Int_t)strip_acc[ilayer].size()) {
                            strip_acc[ilayer][strip_num] += 1.0;
                            if (strip_num < acc_min[ilayer])
                                acc_min[ilayer] = strip_num;
                            if (strip_num > acc_max[ilayer])
                                acc_max[ilayer] = strip_num;
                        } else {
                            // out of the accumulator range: directly, as in the other methods
                            cluster_layers[ilayer].AddStrip(strip_num, 1.0);
                        }
                    }
 
                    // Electron points on the readout (generation data) ---------
                    if (FillGenerationData) {
                        XElectronPos.push_back(x_readout);
                        YElectronPos.push_back(y_readout);
                        ElectronSignal.push_back(1.0);
                    }
                    //----------------------------------------------------------
                }
            }
        }
    }
 
    // Standard method: accumulated signals -> clusters (one AddStrip per strip)
    if (!useANN) {
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            for (Int_t istrip = acc_min[ilayer]; istrip <= acc_max[ilayer]; ++istrip) {
                if (strip_acc[ilayer][istrip] > 0.0)
                    cluster_layers[ilayer].AddStrip(istrip, strip_acc[ilayer][istrip]);
            }
        }
    }
    //----------------------------------------------------------------------
 
    // ANN generation (fdeep) ---------------------------------------------
    if (useANN) {
 
        // Segments for the model: external or from the avalanche of this point
        const vector<Int_t>& seg_used = use_external_segs ? segs : seg_local;
 
        // Normalization (px,py,pz), where r_vec must be equal to 1
        Double_t px_norm = px / r_vec;
        Double_t py_norm = py / r_vec;
        Double_t pz_norm = pz / r_vec;
 
        // Inversion for z-reverse direction particles: the same track segment passed from (x_out, y_out)
        if (pz < 0.0) {
            px_norm = -px_norm;
            py_norm = -py_norm;
            pz_norm = -pz_norm;
        }
 
        Double_t drift_direction = 0;   // 0 is ForwardZAxisEDrift
        if (ElectronDriftDirection == BackwardZAxisEDrift)
            drift_direction = 1;
 
        // Magnetic field: given explicitly or from the current medium configuration
        const Double_t bfield_used = (bfield < 0.0) ? BmnGemStripMedium::GetInstance().BField : bfield;
 
        // Normalize seg. vector: [0, max] -> [0, 1] (max over the whole vector)
        Double_t max_segment_val = *(std::max_element(seg_used.begin(), seg_used.end()));
 
        // Total charge in segment vector
        Double_t total_charge = 0.0;
        for (size_t i = 0; i < seg_used.size(); ++i)
            total_charge += seg_used[i];
 
        // Window center point (x_c, y_c): computed above with GetANNWindowCenterPoint()
 
        // Common part of the model input: [px, py, pz, drift, bfield, s0..s{n_segs_features-1}]
        std::vector<float> input_data;
        input_data.reserve(n_ann_inputs);
        input_data.push_back(static_cast<float>(px_norm));
        input_data.push_back(static_cast<float>(py_norm));
        input_data.push_back(static_cast<float>(pz_norm));
        input_data.push_back(static_cast<float>(drift_direction));
        input_data.push_back(static_cast<float>(bfield_used));
        for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
            input_data.push_back((max_segment_val > 0.0) ? static_cast<float>(seg_used[iseg] / max_segment_val) : 0.0f);
 
        if (total_charge > 0.0) {
            // Strip model: a window of n_ann_strips strips around the window center strip in each layer
            const Int_t center_strip_index = n_ann_strips / 2;   // center strip of the window
 
            // [PARALLEL-ANN] the body of the former loop over layers: one layer, ilayer - its index
            auto process_layer = [&](Int_t ilayer) {
                // The window is built around the center point (middle of the drift gap): it must be inside this layer
                if (!StripLayers[ilayer].IsPointInsideStripLayer(x_c, y_c))
                    return;   // [PARALLEL-ANN] was continue
 
                Int_t n_strips_in_layer = StripLayers[ilayer].GetNStrips();
                if (n_strips_in_layer <= 0)
                    return;   // [PARALLEL-ANN] was continue
 
                // Center point position in strips: integer part - window center strip, fractional part - center_frac
                Double_t center_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_c, y_c);
                Int_t center_strip = (Int_t)TMath::Floor(center_pos);
                Double_t center_frac = center_pos - center_strip;
 
                Double_t strip_layer = static_cast<Double_t>(StripLayers[ilayer].GetType());   // 0 - lower, 1 - upper
                Double_t strip_order =
                    static_cast<Double_t>(StripLayers[ilayer].GetStripNumberingOrder());   // 0 - LeftToRight
 
                // Layer part of the model input: [strip_layer, strip_order, center_frac]
                std::vector<float> layer_input = input_data;
                layer_input.push_back(static_cast<float>(strip_layer));
                layer_input.push_back(static_cast<float>(strip_order));
                layer_input.push_back(static_cast<float>(center_frac));
 
                // Generate strip signals with AI model: n_ann_strips values in [0, 1] (max ~ 1)
                const auto strip_values = ann_model->Predict(layer_input);
                if (strip_values.size() != static_cast<std::size_t>(n_ann_strips)) {
                    LOG(error) << "BmnGemStripModule::"
                               << "AddRealPointFull_ANNStripGenerator_Fdeep"
                               << "(): model output size = " << strip_values.size() << ", expected " << n_ann_strips;
                    return;   // [PARALLEL-ANN] was continue
                }
 
                // Sum of the model output over all strips of the window
                Double_t sum_norm = 0.0;
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip)
                    sum_norm += strip_values.at(istrip);
 
                // Scale coeff for total signal distribution (separate for lower and upper layers)
                const Double_t ann_correction_lower = 0.8;   // Set to 1.0 to disable correction
                const Double_t ann_correction_upper = 0.6;   // Set to 1.0 to disable correction
                const Double_t ann_correction =
                    (StripLayers[ilayer].GetType() == LowerStripLayer) ? ann_correction_lower : ann_correction_upper;
 
                Double_t signal_scale = 0.0;
                if (sum_norm > 0.0)
                    signal_scale = (total_charge / sum_norm) * ann_correction;
 
                Double_t threshold = 0.0001 * total_charge;
 
                // Strips are written directly: window index -> strip number
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip) {
                    Double_t content = strip_values.at(istrip) * signal_scale;
                    if (content <= threshold)
                        continue;
 
                    Int_t strip_num = center_strip - center_strip_index + istrip;
                    if (strip_num < 0 || strip_num >= n_strips_in_layer)
                        continue;   // the window goes out of the layer
 
                    cluster_layers[ilayer].AddStrip(strip_num, content);
                }
            };   // end of process_layer
 
            // [PARALLEL-ANN] layers 1..n-1 in separate threads, layer 0 in the current thread
            // (each layer writes only to its own cluster -> the same result as the sequential loop)
            if (kParallelANNLayers && n_strip_layers > 1) {
                std::vector<std::thread> layer_threads;
                for (Int_t ilayer = 1; ilayer < n_strip_layers; ++ilayer)
                    layer_threads.emplace_back(process_layer, ilayer);
                process_layer(0);
                for (auto& th : layer_threads)
                    th.join();
            } else {
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
                    process_layer(ilayer);
            }
        }
    }
    //----------------------------------------------------------------------
 
    // Segmented charge distribution (generation data) ------------------------
    // (seg_used above is not used after this point)
    if (FillGenerationData) {
        if (use_external_segs)
            SegmentedChargeDistribution = std::move(segs);
        else
            SegmentedChargeDistribution = std::move(seg_local);
    }
    //----------------------------------------------------------------------
 
    // Calculate cluster parameters --------------------------------------------
    vector<Double_t> cluster_mean_position(n_strip_layers, 0.0);
    vector<Double_t> cluster_total_signal(n_strip_layers, 0.0);
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
        if (cluster_layers[ilayer].GetClusterSize() > 0) {
 
            for (int i = 0; i < cluster_layers[ilayer].GetClusterSize(); ++i) {
                cluster_mean_position[ilayer] += (cluster_layers[ilayer].Strips[i] + 0.5)
                                                 * cluster_layers[ilayer].Signals[i];   // as sum of all positions
                cluster_total_signal[ilayer] += cluster_layers[ilayer].Signals[i];
            }
            cluster_mean_position[ilayer] /= cluster_total_signal[ilayer];
        }
 
        cluster_layers[ilayer].MeanPosition = cluster_mean_position[ilayer];   // mean cluster position
        cluster_layers[ilayer].TotalSignal = cluster_total_signal[ilayer];     // total signal of the cluster
    }
 
    if (ElectronDriftDirection == ForwardZAxisEDrift) {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        }
    } else {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        cluster_layers[ilayer].PositionResidual =
            cluster_layers[ilayer].MeanPosition - cluster_layers[ilayer].OriginPosition;
        cluster_layers[ilayer].SetType(StripLayers[ilayer].GetType());
    }
    //----------------------------------------------------------------------
 
    // Added clusters (generation data) ----------------------------------------
    if (FillGenerationData) {
        if (AddedClusters.size() == 0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                AddedClusters.push_back(vector<StripCluster>());
            }
        }
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
        }
    }
    //------------------------------------------------------------------------------
 
    // Add the correct clusters to the strip layers -------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
        }
    }
    //----------------------------------------------------------------------
 
    // Fill strip matches ---------------------------------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                    refID);
        }
    }
    //----------------------------------------------------------------------
 
    RealPointsX.push_back(x);
    RealPointsY.push_back(y);
    RealPointsXout.push_back(x_out);
    RealPointsYout.push_back(y_out);
    RealPointsMC.push_back(refID);
    RealPointsXcenter.push_back(x_c);   // ANN window center point
    RealPointsYcenter.push_back(y_c);
 
    return true;
}
//------------------------------------------------------------------------------
 
// The same as AddRealPointFull_ANNStripGenerator_Fdeep, the standard method is parallel
Bool_t BmnGemStripModule::AddRealPointFull_ANNStripGenerator_Fdeep_Parallel(Double_t x,
                                                                            Double_t y,
                                                                            Double_t z,
                                                                            Double_t px,
                                                                            Double_t py,
                                                                            Double_t pz,
                                                                            Double_t signal,
                                                                            Int_t refID,
                                                                            Int_t pdg_code,
                                                                            Double_t bfield,
                                                                            Bool_t onlyStandard,
                                                                            vector<Int_t> segs)
{
    const Double_t maxThetaAngle = kANNMaxThetaAngle;   // max |theta| (deg) of the track for the ANN
    const Bool_t kParallelANNLayers = true;   // [PARALLEL-ANN] true - ANN prediction for the strip layers in parallel
 
    // Model format: [px, py, pz, drift, bfield, s0..s63, strip_layer, strip_order, center_frac] -> window of N
    // strips (N odd, from the model) centered at the middle of the drift gap (GetANNWindowCenterPoint)
    const Int_t n_track_segments = 64;                    // track segmentation (segmented charge distribution)
    const Int_t n_segs_features = n_track_segments;       // segments passed to the model
    const Int_t n_ann_inputs = 5 + n_segs_features + 3;   // total number of input features
    Int_t n_ann_strips = 0;   // model output: strip window size (odd), taken from the model below
 
    // Test outputs are cleared always (cheap if they are empty)
    ResetElectronPointsAndClusters();
 
    // fdeep model: shared between modules, checked once for each model -------
    std::shared_ptr<const BmnGemStripFdeepModel> ann_holder;
    const BmnGemStripFdeepModel* ann_model = nullptr;
    if (!onlyStandard) {
        ann_holder = BmnGemStripFdeepModel::GetShared();
        ann_model = ann_holder.get();
 
        static std::mutex checked_mutex;
        static std::map<const BmnGemStripFdeepModel*,
                        std::pair<std::shared_ptr<const BmnGemStripFdeepModel>, std::string>>
            checked;   // model -> (holder, why it is not suitable: empty if suitable)
 
        std::lock_guard<std::mutex> lock(checked_mutex);
        if (!ann_model) {
            LOG(warn) << "BmnGemStripModule::" << __func__
                      << "(): fdeep model is not available, the standard method is used";
        } else {
            auto it = checked.find(ann_model);
            if (it == checked.end()) {
                std::string why;
                if (ann_model->GetInputDim() != n_ann_inputs)
                    why = "input " + std::to_string(ann_model->GetInputDim()) + ", expected "
                          + std::to_string(n_ann_inputs);
                else
                    why = CheckANNStripWindowSize(ann_model->GetOutputDim());
                if (why.empty())
                    LOG(info) << "BmnGemStripModule::" << __func__ << "(): fdeep model " << ann_model->GetFilePath()
                              << " is used (ANN for |theta| < " << maxThetaAngle << " deg, window "
                              << ann_model->GetOutputDim() << " strips)";
                it = checked.emplace(ann_model, std::make_pair(ann_holder, why)).first;
            }
            if (!it->second.second.empty()) {
                LOG(warn) << "BmnGemStripModule::" << __func__ << "(): the model " << ann_model->GetFilePath()
                          << " is not suitable (" << it->second.second << "), the standard method is used";
                ann_model = nullptr;
            }
        }
    }
    if (ann_model)
        n_ann_strips = ann_model->GetOutputDim();   // window size of the model (checked above: odd, >= 3)
    //--------------------------------------------------------------------------
 
    // Condition: a start point is inside the module and outside its dead zone (in strip layers)
    if (!IsPointInsideModule(x, y)) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
 
    // Particle: pdg_code = 0 -> from the MC track of the point refID (as in AddRealPointFull) --
    Int_t pdg = pdg_code;
    if (pdg == 0) {
        FairRunSim* run = FairRunSim::Instance();
        BmnGemStripDigitizer* digi = run ? (BmnGemStripDigitizer*)run->GetTask("GemDigitizer") : nullptr;
        TClonesArray* points = digi ? digi->GetMCPoints() : nullptr;
        TClonesArray* mctracks = digi ? digi->GetMCTracks() : nullptr;
 
        CbmMCTrack* mcTr = nullptr;
        if (points && mctracks && refID >= 0 && refID < points->GetEntriesFast()) {
            CbmStsPoint* point = (CbmStsPoint*)points->UncheckedAt(refID);
            if (point && point->GetTrackID() >= 0 && point->GetTrackID() < mctracks->GetEntriesFast())
                mcTr = (CbmMCTrack*)mctracks->UncheckedAt(point->GetTrackID());
        }
 
        if (mcTr) {
            pdg = mcTr->GetPdgCode();
        } else {
            static std::atomic<bool> no_mc_reported{false};
            if (!no_mc_reported.exchange(true))
                LOG(warn) << "BmnGemStripModule::" << __func__
                          << "(): pdg_code = 0 and the MC track is not available (no GemDigitizer or MC point), "
                             "the He4 mass is used (this warning is printed once)";
        }
    }
 
    TVector3 mom(px, py, pz);
    TParticlePDG* part = (pdg != 0) ? TDatabasePDG::Instance()->GetParticle(pdg) : nullptr;
 
    Double_t mnoc = 1.0, mass = 3.7283999;   // He4
    if (part)
        mass = part->Mass();
    if (mass > 0.0001) {
        TLorentzVector lorv;
        lorv.SetVectM(mom, mass);
        Double_t beta = lorv.Beta();
        Double_t gamma = lorv.Gamma();
        Double_t charge = (part) ? part->Charge() / 3 : 1;
        mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
        if (mnoc < 0.5)
            return kFALSE;
        if (MaxClustersPerCm > 0.0 && mnoc > MaxClustersPerCm)
            mnoc = MaxClustersPerCm;
    }
    //--------------------------------------------------------------------------
 
    // Distance from entry point (x,y,z) to exit point (along z-axis)
    Double_t z_distance_from_in_to_out;
 
    // Condition: track direction along z-axis
    if (pz > 0.0) {
        z_distance_from_in_to_out = ModuleThickness - (z - ZStartModulePosition);
    } else {
        z_distance_from_in_to_out = z - ZStartModulePosition;
    }
 
    // Condition: track distance along z-axis must be not zero
    if (z_distance_from_in_to_out <= 0.0) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") has a track with an incorrect length";
        return false;
    }
 
    // Scale coefficient of the track vector (the path from entry point to exit point)
    Double_t vector_coeff = fabs(z_distance_from_in_to_out / pz);
 
    // Components of the track vector
    Double_t x_vec_from_in_to_out = vector_coeff * px;
    Double_t y_vec_from_in_to_out = vector_coeff * py;
    Double_t z_vec_from_in_to_out = vector_coeff * pz;
 
    // Exit point coordinates (x_out, y_out, z_out)
    Double_t x_out = x + x_vec_from_in_to_out;
    Double_t y_out = y + y_vec_from_in_to_out;
    Double_t z_out = z + z_vec_from_in_to_out;
 
    // Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
    // Condition: x-coordinate of the exit point is inside the module
    if (x_out < XMinModule || x_out > XMaxModule) {
        if (x_out < XMinModule)
            x_out = XMinModule;
        if (x_out > XMaxModule)
            x_out = XMaxModule;
 
        x_vec_from_in_to_out = x_out - x;
        vector_coeff = x_vec_from_in_to_out / px;
 
        y_vec_from_in_to_out = vector_coeff * py;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        y_out = y + y_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    // Condition: y-coordinate of the exit point is inside the module
    if (y_out < YMinModule || y_out > YMaxModule) {
        if (y_out < YMinModule)
            y_out = YMinModule;
        if (y_out > YMaxModule)
            y_out = YMaxModule;
 
        y_vec_from_in_to_out = y_out - y;
        vector_coeff = y_vec_from_in_to_out / py;
 
        x_vec_from_in_to_out = vector_coeff * px;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        x_out = x + x_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    //--------------------------------------------------------------------------
 
    // Center point of the ANN strip window: the track at the middle of the drift gap (computed always, also with
    // onlyStandard: it is stored with the point, see GetRealPointXcenter(), for the training samples)
    Double_t x_c = x, y_c = y;
    GetANNWindowCenterPoint(x, y, z, x_out, y_out, z_out, x_c, y_c);
 
    // Common track length
    Double_t track_length =
        std::sqrt((x_vec_from_in_to_out) * (x_vec_from_in_to_out) + (y_vec_from_in_to_out) * (y_vec_from_in_to_out)
                  + (z_vec_from_in_to_out) * (z_vec_from_in_to_out));
 
    // Which method is used for this point: standard (electrons -> strips) or ANN (model -> strips)
    // (without a suitable model or with onlyStandard the standard method is used)
    Double_t r_vec = TMath::Sqrt(px * px + py * py + pz * pz);
    Double_t theta_deg_abs = TMath::ACos(TMath::Abs(pz) / r_vec) * TMath::RadToDeg();
    const Bool_t useANN = (theta_deg_abs < maxThetaAngle) && ann_model;
 
    // ANN with external segments: the avalanche is not needed
    const Bool_t use_external_segs = useANN && !segs.empty();
    if (use_external_segs && segs.size() < static_cast<size_t>(n_segs_features)) {
        LOG(error) << "BmnGemStripModule::" << __func__ << "(): segment vector size = " << segs.size()
                   << ", the model requires " << n_segs_features;
        return false;
    }
 
    // Each level - distance to the readout plane
    Double_t level1 = InductionGapThickness;
    Double_t level2 = InductionGapThickness + SecondTransferGapThickness;
    Double_t level3 = InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness;
 
    // Mean electron shift along x-axis (under the influence of the Lorentz force)
    Double_t xmean;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
 
    // Sigma electron smearing
    Double_t sigma;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3+ p4*x^4 + p5*x^5)
 
    // Reading the medium configuration -------------------------------------
    Double_t p0_xmean = BmnGemStripMedium::GetInstance().p0_xmean;
    Double_t p1_xmean = BmnGemStripMedium::GetInstance().p1_xmean;
    Double_t p2_xmean = BmnGemStripMedium::GetInstance().p2_xmean;
    Double_t p3_xmean = BmnGemStripMedium::GetInstance().p3_xmean;
 
    Double_t p0_sigma = BmnGemStripMedium::GetInstance().p0_sigma;
    Double_t p1_sigma = BmnGemStripMedium::GetInstance().p1_sigma;
    Double_t p2_sigma = BmnGemStripMedium::GetInstance().p2_sigma;
    Double_t p3_sigma = BmnGemStripMedium::GetInstance().p3_sigma;
    Double_t p4_sigma = BmnGemStripMedium::GetInstance().p4_sigma;
    Double_t p5_sigma = BmnGemStripMedium::GetInstance().p5_sigma;
    //----------------------------------------------------------------------
 
    // Initialization of a strip cluster in each strip layer
    Int_t n_strip_layers = StripLayers.size();   // number of layers in the module
    vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
    // Segmented charge distribution of this point (always computed when the avalanche is simulated)
    vector<Int_t> seg_local(n_track_segments, 0);
 
    // [PARALLEL-STD] Standard method: parallel processing of the avalanche electrons -----
    // true  - electron (x,y) are taken from gRandom sequentially, as in AddRealPointFull_ANNStripGenerator_Fdeep:
    //         bitwise identical result, only the strip assignment is parallel;
    // false - electron (x,y) are generated in parallel by per-task generators (seeded from gRandom):
    //         faster, statistically equivalent to the sequential methods, bitwise identical to
    //         AddRealPointFull_ANNStripGenerator_MLP_Parallel (the same scheme of the random numbers).
    constexpr Bool_t kExactGRandomSequence = false;
    constexpr Long64_t kElectronsPerTask = 4096;           // electrons per parallel task
    constexpr Long64_t kMinElectronsForParallel = 20000;   // below this - sequential execution
 
    struct AvalancheTask   // used only if !kExactGRandomSequence
    {
        Double_t xcenter;   // mean x on the readout (with the Lorentz shift)
        Double_t ycenter;   // mean y on the readout
        Double_t sigma;     // smearing
        Long64_t offset;    // index of the first electron of the task (from first_electron)
        Long64_t n;         // number of electrons in the task
    };
    std::vector<AvalancheTask> aval_tasks;
    Long64_t n_aval_electrons = 0;
    std::vector<Double_t> exact_x, exact_y;   // electron positions, used only if kExactGRandomSequence
 
    // Persistent thread pool: the threads are created once (on the first call) and wait for work,
    // so no threads are created/joined in each call. Lives until the end of the program.
    // Work items are taken dynamically (atomic counter) -> load balancing for avalanches of different size.
    struct ThreadPool
    {
        std::vector<std::thread> workers;
        std::mutex mtx;                     // protects the job description and the counters below
        std::mutex run_mtx;                 // one job at a time (in case of calls from several threads)
        std::condition_variable cv_start;   // workers: a new job is available
        std::condition_variable cv_done;    // caller: all workers have finished the job
        const std::function<void(size_t)>* job = nullptr;
        size_t n_items = 0;
        std::atomic<size_t> next{0};   // next work item
        size_t generation = 0;         // job counter
        size_t n_busy = 0;             // workers that have not finished the current job yet
        bool stop = false;
 
        explicit ThreadPool(size_t n_workers)
        {
            workers.reserve(n_workers);
            for (size_t i = 0; i < n_workers; ++i)
                workers.emplace_back([this] { Loop(); });
        }
 
        ~ThreadPool()
        {
            {
                std::lock_guard<std::mutex> lk(mtx);
                stop = true;
            }
            cv_start.notify_all();
            for (auto& w : workers)
                w.join();
        }
 
        void Work()
        {
            for (size_t i = next.fetch_add(1); i < n_items; i = next.fetch_add(1))
                (*job)(i);
        }
 
        void Loop()
        {
            size_t seen = 0;
            for (;;) {
                {
                    std::unique_lock<std::mutex> lk(mtx);
                    cv_start.wait(lk, [&] { return stop || generation != seen; });
                    if (stop)
                        return;
                    seen = generation;
                }
                Work();
                {
                    std::lock_guard<std::mutex> lk(mtx);
                    if (--n_busy == 0)
                        cv_done.notify_one();
                }
            }
        }
 
        // Runs func(i) for i in [0, n) on the pool threads and the calling thread
        void Run(size_t n, const std::function<void(size_t)>& func)
        {
            if (n == 0)
                return;
            if (workers.empty() || n == 1) {
                for (size_t i = 0; i < n; ++i)
                    func(i);
                return;
            }
 
            std::lock_guard<std::mutex> run_lock(run_mtx);
            {
                std::lock_guard<std::mutex> lk(mtx);
                job = &func;
                n_items = n;
                next.store(0);
                n_busy = workers.size();
                ++generation;
            }
            cv_start.notify_all();
 
            Work();   // the calling thread works too
 
            std::unique_lock<std::mutex> lk(mtx);
            cv_done.wait(lk, [&] { return n_busy == 0; });
            job = nullptr;
        }
    };
 
    // hardware threads - 1 workers (the calling thread is the last one)
    static ThreadPool pool(std::max(1u, std::thread::hardware_concurrency()) - 1);
 
    auto parallel_for = [](size_t n, const std::function<void(size_t)>& func) { pool.Run(n, func); };
    //--------------------------------------------------------------------------
 
    if (!use_external_segs) {
 
        Double_t current_length = 0.0;         // traversed length (counter)
        Double_t current_length_ratio = 0.0;   // ratio of the traversed length to common track length (counter)
        Double_t current_step = 0.0;           // current distance between two collision points
 
        // Collection of collision points
        std::vector<CollPoint> collision_points;
 
        // Reading the MCD value (mean collision distance = mean free flight path) [cm]
        MCD = BmnGemStripMedium::GetInstance().MCD;
        MCD = 1 / mnoc;   // AZ-110124
 
        while (current_length < track_length) {
 
            current_step = gRandom->Exp(MCD);
            current_length += current_step;
 
            if (current_length > track_length)
                break;
 
            current_length_ratio = current_length / track_length;
 
            Double_t current_x = x + current_length_ratio * x_vec_from_in_to_out;
            Double_t current_y = y + current_length_ratio * y_vec_from_in_to_out;
            Double_t current_z = z + current_length_ratio * z_vec_from_in_to_out;
 
            collision_points.push_back(CollPoint(current_x, current_y, current_z));
        }
 
        // Electron avalanche producing for each collision point at the track
        // (in the ANN mode only the segmented charge distribution is filled here)
        for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
 
            Double_t xcoll = collision_points[icoll].x;
            Double_t ycoll = collision_points[icoll].y;
            Double_t zcoll = collision_points[icoll].z;
 
            Double_t zdist;   // current z-distance to the readout
 
            // Find z-distance to the readout depending on the electron drift direction
            if (ElectronDriftDirection == ForwardZAxisEDrift) {
                zdist = (ModuleThickness + ZStartModulePosition) - zcoll;
            } else {
                zdist = zcoll - ZStartModulePosition;
            }
            //------------------------------------------------------------------
 
            // polynomials
            xmean = p0_xmean + p1_xmean * zdist + p2_xmean * zdist * zdist + p3_xmean * zdist * zdist * zdist;
            sigma = p0_sigma + p1_sigma * zdist + p2_sigma * zdist * zdist + p3_sigma * zdist * zdist * zdist
                    + p4_sigma * zdist * zdist * zdist * zdist + p5_sigma * zdist * zdist * zdist * zdist * zdist;
 
            // Number of electrons in the current collision
            Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
            if (n_electrons_cluster < 1)
                n_electrons_cluster = 1;   // min
            if (n_electrons_cluster > 6)
                n_electrons_cluster = 6;   // max
 
            // Filling charge distribution -------------------------------------
            Double_t dist_from_start = 0.0;
            if (pz > 0.0) {
                dist_from_start =
                    std::sqrt((xcoll - x) * (xcoll - x) + (ycoll - y) * (ycoll - y) + (zcoll - z) * (zcoll - z));
            } else {
                dist_from_start = std::sqrt((xcoll - x_out) * (xcoll - x_out) + (ycoll - y_out) * (ycoll - y_out)
                                            + (zcoll - z_out) * (zcoll - z_out));
            }
            Int_t segment = (Int_t)(dist_from_start / track_length * n_track_segments);
            if (segment >= n_track_segments)
                segment = n_track_segments - 1;
            //------------------------------------------------------------------
 
            for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
 
                // Electron gain in each GEM cascade (three draws always, as in the other methods)
                Double_t gain_gem1 = gRandom->Exp(15);
                Double_t gain_gem2 = gRandom->Exp(15);
                Double_t gain_gem3 = gRandom->Exp(15);
 
                if (gain_gem1 < 1.0)
                    gain_gem1 = 1.0;
                if (gain_gem2 < 1.0)
                    gain_gem2 = 1.0;
                if (gain_gem3 < 1.0)
                    gain_gem3 = 1.0;
 
                int total_gain = 0;
 
                if (zdist < level1) {
                    total_gain = 1.0;
                }
                if (zdist >= level1 && zdist < level2) {
                    total_gain = gain_gem3;
                }
                if (zdist >= level2 && zdist < level3) {
                    total_gain = gain_gem3 * gain_gem2;
                }
                if (zdist >= level3) {
                    total_gain = gain_gem3 * gain_gem2 * gain_gem1;
                }
 
                seg_local[segment] += total_gain;
 
                // ANN mode: electrons are not projected on the readout, strips are taken from the model
                if (useANN)
                    continue;
 
                // Standard method ---------------------------------------------
                // [PARALLEL-STD] Only the electron positions (exact mode) or the tasks (fast mode) are produced
                // here; strips are assigned later in parallel (see "Standard method: parallel strip assignment")
                // x-shift of the electon depends on the electron drift direction
                const Double_t x_center =
                    (ElectronDriftDirection == ForwardZAxisEDrift) ? (xcoll + xmean) : (xcoll - xmean);
 
                if (kExactGRandomSequence) {
                    // the same gRandom calls in the same order as in AddRealPointFull_ANNStripGenerator_Fdeep
                    for (int igain = 0; igain < total_gain; ++igain) {
                        exact_x.push_back(gRandom->Gaus(x_center, sigma));
                        exact_y.push_back(gRandom->Gaus(ycoll, sigma));
                    }
                } else {
                    // positions are generated later in parallel, here - only the tasks
                    for (Long64_t done = 0; done < total_gain; done += kElectronsPerTask) {
                        aval_tasks.push_back({x_center, ycoll, sigma, n_aval_electrons + done,
                                              std::min<Long64_t>(kElectronsPerTask, total_gain - done)});
                    }
                    n_aval_electrons += total_gain;
                }
            }
        }
    }
 
    // [PARALLEL-STD] Standard method: parallel strip assignment ------------
    // Each task counts its electrons per strip (compact local counters, no per-electron storage).
    // The counters of the tasks are summed in the task order and added to the clusters with one AddStrip()
    // per strip (as in AddRealPointFull_ANNStripGenerator_MLP): the signals are integer numbers of electrons,
    // so the clusters are bitwise the same as with one AddStrip() per electron and do not depend on the threads.
    if (!useANN) {
 
        // Number of avalanche electrons of this point
        const size_t n_electrons = kExactGRandomSequence ? exact_x.size() : static_cast<size_t>(n_aval_electrons);
 
        if (n_electrons > 0) {
 
            const Bool_t run_parallel = static_cast<Long64_t>(n_electrons) >= kMinElectronsForParallel;
            const size_t chunk_size = static_cast<size_t>(kElectronsPerTask);
 
            // Layer invariant (computed once instead of for each electron)
            std::vector<char> layer_vertical(n_strip_layers);
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
                layer_vertical[ilayer] = TMath::Abs(StripLayers[ilayer].GetAngleDeg()) < 1.0;
 
            // Counters of one task: in each layer the strips [lo, lo + cnt.size()) (the range grows if needed)
            struct TaskCounts
            {
                std::vector<Int_t> lo;
                std::vector<std::vector<Int_t>> cnt;
            };
            const size_t n_tasks =
                kExactGRandomSequence ? (n_electrons + chunk_size - 1) / chunk_size : aval_tasks.size();
            std::vector<TaskCounts> task_counts(n_tasks);
 
            // One electron -> +1 on its strip in each layer (the same conditions as in the old code)
            auto count_electron = [&](TaskCounts& tc, Double_t x_readout, Double_t y_readout) {
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    // NOTE: these BmnGemStripLayer methods are called from several threads,
                    //       they must not modify the layer
                    // Condition: end electron position must be inside the layer
                    if (!StripLayers[ilayer].IsPointInsideStripLayer(x_readout, y_readout))
                        continue;
 
                    Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_readout, y_readout);
                    // AZ-100124 - account for strip width
                    Double_t dstrip = strip_pos - (int)strip_pos;
                    if (layer_vertical[ilayer]) {
                        // Vertical strips
                        if (dstrip >= 0.4 && dstrip < 0.6)
                            continue;
                    } else {
                        // Rotated strips
                        if (dstrip >= 0.3 && dstrip < 0.7)
                            continue;
                    }
                    // AZ
 
                    const Int_t strip = (Int_t)strip_pos;
                    std::vector<Int_t>& cnt = tc.cnt[ilayer];
                    Int_t& lo = tc.lo[ilayer];
                    if (cnt.empty()) {
                        lo = strip;
                        cnt.assign(1, 0);
                    } else if (strip < lo) {
                        cnt.insert(cnt.begin(), lo - strip, 0);
                        lo = strip;
                    } else if (strip >= lo + static_cast<Int_t>(cnt.size())) {
                        cnt.resize(strip - lo + 1, 0);
                    }
                    ++cnt[strip - lo];
                }
            };
 
            auto init_counts = [&](TaskCounts& tc) {
                tc.lo.assign(n_strip_layers, 0);
                tc.cnt.assign(n_strip_layers, std::vector<Int_t>());
            };
 
            if (kExactGRandomSequence) {
                // (a) Exact mode: positions are already generated (gRandom, sequentially) -> only geometry, by chunks
                auto process_chunk = [&](size_t ichunk) {
                    TaskCounts& tc = task_counts[ichunk];
                    init_counts(tc);
                    const size_t begin = ichunk * chunk_size;
                    const size_t end = std::min(n_electrons, begin + chunk_size);
                    for (size_t ie = begin; ie < end; ++ie)
                        count_electron(tc, exact_x[ie], exact_y[ie]);
                };
 
                if (run_parallel)
                    parallel_for(n_tasks, process_chunk);
                else
                    for (size_t ichunk = 0; ichunk < n_tasks; ++ichunk)
                        process_chunk(ichunk);
 
                // Electron points on the readout (generation data)
                if (FillGenerationData) {
                    XElectronPos.insert(XElectronPos.end(), exact_x.begin(), exact_x.end());
                    YElectronPos.insert(YElectronPos.end(), exact_y.begin(), exact_y.end());
                    ElectronSignal.resize(ElectronSignal.size() + n_electrons, 1.0);
                }
            } else {
                // (a) Fast mode: generation + geometry in one pass (one thread launch, one pass over memory)
 
                // base seed from gRandom -> reproducible with the seed of gRandom, independent of the threads
                const UInt_t base_seed = static_cast<UInt_t>(gRandom->Integer(0xFFFFFFFFu));
 
                // Electron points on the readout (generation data): each task writes its own range
                const size_t first_electron = XElectronPos.size();   // 0 after ResetElectronPointsAndClusters()
                if (FillGenerationData) {
                    XElectronPos.resize(first_electron + n_electrons);
                    YElectronPos.resize(first_electron + n_electrons);
                    ElectronSignal.resize(first_electron + n_electrons, 1.0);
                }
 
                auto process_task = [&](size_t itask) {
                    const AvalancheTask& t = aval_tasks[itask];
                    TaskCounts& tc = task_counts[itask];
                    init_counts(tc);
                    std::seed_seq seq{base_seed, static_cast<UInt_t>(itask)};
                    std::mt19937_64 rng(seq);
                    std::normal_distribution<Double_t> gaus(0.0, 1.0);   // mean + sigma*N(0,1), as TRandom::Gaus
 
                    for (Long64_t i = 0; i < t.n; ++i) {
                        const Double_t x_readout = t.xcenter + t.sigma * gaus(rng);
                        const Double_t y_readout = t.ycenter + t.sigma * gaus(rng);
 
                        if (FillGenerationData) {
                            const size_t ie = first_electron + static_cast<size_t>(t.offset + i);
                            XElectronPos[ie] = x_readout;
                            YElectronPos[ie] = y_readout;
                        }
 
                        count_electron(tc, x_readout, y_readout);
                    }
                };
 
                if (run_parallel)
                    parallel_for(n_tasks, process_task);
                else
                    for (size_t itask = 0; itask < n_tasks; ++itask)
                        process_task(itask);
            }
 
            // (b) Filling the clusters: the counters of all tasks are summed (in the task order),
            //     one AddStrip() per strip; layers are independent -> in parallel
            auto fill_layer = [&](size_t ilayer) {
                StripCluster& cluster = cluster_layers[ilayer];
                const Int_t n_strips_in_layer = std::max(0, StripLayers[ilayer].GetNStrips());
                std::vector<Double_t> strip_acc(n_strips_in_layer, 0.0);
                Int_t acc_min = std::numeric_limits<Int_t>::max();
                Int_t acc_max = -1;
 
                for (size_t itask = 0; itask < n_tasks; ++itask) {
                    const Int_t lo = task_counts[itask].lo[ilayer];
                    const std::vector<Int_t>& cnt = task_counts[itask].cnt[ilayer];
                    for (size_t i = 0; i < cnt.size(); ++i) {
                        if (cnt[i] == 0)
                            continue;
                        const Int_t strip = lo + static_cast<Int_t>(i);
                        if (strip >= 0 && strip < n_strips_in_layer) {
                            strip_acc[strip] += cnt[i];
                            if (strip < acc_min)
                                acc_min = strip;
                            if (strip > acc_max)
                                acc_max = strip;
                        } else {
                            // out of the accumulator range: directly (the same sum as cnt[i] x AddStrip(strip, 1.0))
                            cluster.AddStrip(strip, static_cast<Double_t>(cnt[i]));
                        }
                    }
                }
 
                for (Int_t istrip = acc_min; istrip <= acc_max; ++istrip) {
                    if (strip_acc[istrip] > 0.0)
                        cluster.AddStrip(istrip, strip_acc[istrip]);
                }
            };
 
            if (run_parallel)
                parallel_for(static_cast<size_t>(n_strip_layers), fill_layer);
            else
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
                    fill_layer(ilayer);
        }
    }
    //----------------------------------------------------------------------
 
    // ANN generation (fdeep) ---------------------------------------------
    if (useANN) {
 
        // Segments for the model: external or from the avalanche of this point
        const vector<Int_t>& seg_used = use_external_segs ? segs : seg_local;
 
        // Normalization (px,py,pz), where r_vec must be equal to 1
        Double_t px_norm = px / r_vec;
        Double_t py_norm = py / r_vec;
        Double_t pz_norm = pz / r_vec;
 
        // Inversion for z-reverse direction particles: the same track segment passed from (x_out, y_out)
        if (pz < 0.0) {
            px_norm = -px_norm;
            py_norm = -py_norm;
            pz_norm = -pz_norm;
        }
 
        Double_t drift_direction = 0;   // 0 is ForwardZAxisEDrift
        if (ElectronDriftDirection == BackwardZAxisEDrift)
            drift_direction = 1;
 
        // Magnetic field: given explicitly or from the current medium configuration
        const Double_t bfield_used = (bfield < 0.0) ? BmnGemStripMedium::GetInstance().BField : bfield;
 
        // Normalize seg. vector: [0, max] -> [0, 1] (max over the whole vector)
        Double_t max_segment_val = *(std::max_element(seg_used.begin(), seg_used.end()));
 
        // Total charge in segment vector
        Double_t total_charge = 0.0;
        for (size_t i = 0; i < seg_used.size(); ++i)
            total_charge += seg_used[i];
 
        // Window center point (x_c, y_c): computed above with GetANNWindowCenterPoint()
 
        // Common part of the model input: [px, py, pz, drift, bfield, s0..s{n_segs_features-1}]
        std::vector<float> input_data;
        input_data.reserve(n_ann_inputs);
        input_data.push_back(static_cast<float>(px_norm));
        input_data.push_back(static_cast<float>(py_norm));
        input_data.push_back(static_cast<float>(pz_norm));
        input_data.push_back(static_cast<float>(drift_direction));
        input_data.push_back(static_cast<float>(bfield_used));
        for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
            input_data.push_back((max_segment_val > 0.0) ? static_cast<float>(seg_used[iseg] / max_segment_val) : 0.0f);
 
        if (total_charge > 0.0) {
            // Strip model: a window of n_ann_strips strips around the window center strip in each layer
            const Int_t center_strip_index = n_ann_strips / 2;   // center strip of the window
 
            // [PARALLEL-STD] the body of the former loop over layers: one layer, ilayer - its index
            auto process_layer = [&](size_t ilayer) {
                // The window is built around the center point (middle of the drift gap): it must be inside this layer
                if (!StripLayers[ilayer].IsPointInsideStripLayer(x_c, y_c))
                    return;   // [PARALLEL-STD] was continue
 
                Int_t n_strips_in_layer = StripLayers[ilayer].GetNStrips();
                if (n_strips_in_layer <= 0)
                    return;   // [PARALLEL-STD] was continue
 
                // Center point position in strips: integer part - window center strip, fractional part - center_frac
                Double_t center_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_c, y_c);
                Int_t center_strip = (Int_t)TMath::Floor(center_pos);
                Double_t center_frac = center_pos - center_strip;
 
                Double_t strip_layer = static_cast<Double_t>(StripLayers[ilayer].GetType());   // 0 - lower, 1 - upper
                Double_t strip_order =
                    static_cast<Double_t>(StripLayers[ilayer].GetStripNumberingOrder());   // 0 - LeftToRight
 
                // Layer part of the model input: [strip_layer, strip_order, center_frac]
                std::vector<float> layer_input = input_data;
                layer_input.push_back(static_cast<float>(strip_layer));
                layer_input.push_back(static_cast<float>(strip_order));
                layer_input.push_back(static_cast<float>(center_frac));
 
                // Generate strip signals with AI model: n_ann_strips values in [0, 1] (max ~ 1)
                const auto strip_values = ann_model->Predict(layer_input);
                if (strip_values.size() != static_cast<std::size_t>(n_ann_strips)) {
                    LOG(error) << "BmnGemStripModule::"
                               << "AddRealPointFull_ANNStripGenerator_Fdeep_Parallel"
                               << "(): model output size = " << strip_values.size() << ", expected " << n_ann_strips;
                    return;   // [PARALLEL-STD] was continue
                }
 
                // Sum of the model output over all strips of the window
                Double_t sum_norm = 0.0;
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip)
                    sum_norm += strip_values.at(istrip);
 
                // Scale coeff for total signal distribution (separate for lower and upper layers)
                const Double_t ann_correction_lower = 0.8;   // Set to 1.0 to disable correction
                const Double_t ann_correction_upper = 0.6;   // Set to 1.0 to disable correction
                const Double_t ann_correction =
                    (StripLayers[ilayer].GetType() == LowerStripLayer) ? ann_correction_lower : ann_correction_upper;
 
                Double_t signal_scale = 0.0;
                if (sum_norm > 0.0)
                    signal_scale = (total_charge / sum_norm) * ann_correction;
 
                Double_t threshold = 0.0001 * total_charge;
 
                // Strips are written directly: window index -> strip number
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip) {
                    Double_t content = strip_values.at(istrip) * signal_scale;
                    if (content <= threshold)
                        continue;
 
                    Int_t strip_num = center_strip - center_strip_index + istrip;
                    if (strip_num < 0 || strip_num >= n_strips_in_layer)
                        continue;   // the window goes out of the layer
 
                    cluster_layers[ilayer].AddStrip(strip_num, content);
                }
            };
 
            // [PARALLEL-STD] layers in parallel (each layer writes only to its own cluster) or sequentially
            if (kParallelANNLayers)
                parallel_for(static_cast<size_t>(n_strip_layers), process_layer);
            else
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
                    process_layer(ilayer);
        }
    }
    //----------------------------------------------------------------------
 
    // Segmented charge distribution (generation data) ------------------------
    // (seg_used above is not used after this point)
    if (FillGenerationData) {
        if (use_external_segs)
            SegmentedChargeDistribution = std::move(segs);
        else
            SegmentedChargeDistribution = std::move(seg_local);
    }
    //----------------------------------------------------------------------
 
    // Calculate cluster parameters --------------------------------------------
    vector<Double_t> cluster_mean_position(n_strip_layers, 0.0);
    vector<Double_t> cluster_total_signal(n_strip_layers, 0.0);
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
        if (cluster_layers[ilayer].GetClusterSize() > 0) {
 
            for (int i = 0; i < cluster_layers[ilayer].GetClusterSize(); ++i) {
                cluster_mean_position[ilayer] += (cluster_layers[ilayer].Strips[i] + 0.5)
                                                 * cluster_layers[ilayer].Signals[i];   // as sum of all positions
                cluster_total_signal[ilayer] += cluster_layers[ilayer].Signals[i];
            }
            cluster_mean_position[ilayer] /= cluster_total_signal[ilayer];
        }
 
        cluster_layers[ilayer].MeanPosition = cluster_mean_position[ilayer];   // mean cluster position
        cluster_layers[ilayer].TotalSignal = cluster_total_signal[ilayer];     // total signal of the cluster
    }
 
    if (ElectronDriftDirection == ForwardZAxisEDrift) {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        }
    } else {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        cluster_layers[ilayer].PositionResidual =
            cluster_layers[ilayer].MeanPosition - cluster_layers[ilayer].OriginPosition;
        cluster_layers[ilayer].SetType(StripLayers[ilayer].GetType());
    }
    //----------------------------------------------------------------------
 
    // Added clusters (generation data) ----------------------------------------
    if (FillGenerationData) {
        if (AddedClusters.size() == 0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                AddedClusters.push_back(vector<StripCluster>());
            }
        }
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
        }
    }
    //------------------------------------------------------------------------------
 
    // Add the correct clusters to the strip layers -------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
        }
    }
    //----------------------------------------------------------------------
 
    // Fill strip matches ---------------------------------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                    refID);
        }
    }
    //----------------------------------------------------------------------
 
    RealPointsX.push_back(x);
    RealPointsY.push_back(y);
    RealPointsXout.push_back(x_out);
    RealPointsYout.push_back(y_out);
    RealPointsMC.push_back(refID);
    RealPointsXcenter.push_back(x_c);   // ANN window center point
    RealPointsYcenter.push_back(y_c);
 
    return true;
}
//------------------------------------------------------------------------------
 
// Strip generation with the own MLP inference (BmnGemStripMLPModel): standard method or ANN (|theta| <
// kANNMaxThetaAngle). pdg_code = 0 -> from MC, bfield < 0 -> from BmnGemStripMedium; ANN with external segs - no
// avalanche. Generation data only with SetFillGenerationData(kTRUE)
Bool_t BmnGemStripModule::AddRealPointFull_ANNStripGenerator_MLP(Double_t x,
                                                                 Double_t y,
                                                                 Double_t z,
                                                                 Double_t px,
                                                                 Double_t py,
                                                                 Double_t pz,
                                                                 Double_t signal,
                                                                 Int_t refID,
                                                                 Int_t pdg_code,
                                                                 Double_t bfield,
                                                                 Bool_t onlyStandard,
                                                                 vector<Int_t> segs)
{
    const Double_t maxThetaAngle = kANNMaxThetaAngle;   // max |theta| (deg) of the track for the ANN
 
    // Model format: [px, py, pz, drift, bfield, s0..s63, strip_layer, strip_order, center_frac] -> window of N
    // strips (N odd, from the model) centered at the middle of the drift gap (GetANNWindowCenterPoint)
    const Int_t n_track_segments = 64;                   // track segmentation (segmented charge distribution)
    const Int_t n_segs_features = n_track_segments;      // segments passed to the model
    const Int_t n_common_inputs = 5 + n_segs_features;   // the same for all strip layers of the point
    const Int_t n_layer_inputs = 3;                      // [strip_layer, strip_order, center_frac]
    Int_t n_ann_strips = 0;   // model output: strip window size (odd), taken from the model below
 
    // Test outputs are cleared always (cheap if they are empty)
    ResetElectronPointsAndClusters();
 
    // MLP model: shared between modules, checked once for each model ---------
    std::shared_ptr<const BmnGemStripMLPModel> mlp_holder;
    const BmnGemStripMLPModel* mlp_model = nullptr;
    if (!onlyStandard) {
        mlp_holder = BmnGemStripMLPModel::GetShared();
        mlp_model = mlp_holder.get();
 
        static std::mutex checked_mutex;
        static std::map<const BmnGemStripMLPModel*, std::pair<std::shared_ptr<const BmnGemStripMLPModel>, std::string>>
            checked;   // model -> (holder, why it is not suitable: empty if suitable)
 
        std::lock_guard<std::mutex> lock(checked_mutex);
        if (!mlp_model) {
            LOG(warn) << "BmnGemStripModule::" << __func__
                      << "(): MLP model is not available, the standard method is used";
        } else {
            auto it = checked.find(mlp_model);
            if (it == checked.end()) {
                std::vector<std::string> expected = {"px", "py", "pz", "drift", "bfield"};
                for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
                    expected.push_back("s" + std::to_string(iseg));
                expected.insert(expected.end(), {"strip_layer", "strip_order", "center_frac"});
 
                std::string why;
                Bool_t ok = mlp_model->CheckFeatureNames(expected, &why);
                if (ok) {
                    why = CheckANNStripWindowSize(mlp_model->GetOutputDim());
                    ok = why.empty();
                }
                if (ok) {
                    why.clear();
                    LOG(info) << "BmnGemStripModule::" << __func__ << "(): MLP model " << mlp_model->GetFilePath()
                              << " is used (ANN for |theta| < " << maxThetaAngle << " deg, window "
                              << mlp_model->GetOutputDim() << " strips)";
                } else if (why.empty())
                    why = "unknown reason";
                it = checked.emplace(mlp_model, std::make_pair(mlp_holder, why)).first;
            }
            if (!it->second.second.empty()) {
                LOG(warn) << "BmnGemStripModule::" << __func__ << "(): the model " << mlp_model->GetFilePath()
                          << " is not suitable (" << it->second.second << "), the standard method is used";
                mlp_model = nullptr;
            }
        }
    }
    if (mlp_model)
        n_ann_strips = mlp_model->GetOutputDim();   // window size of the model (checked above: odd, >= 3)
    //--------------------------------------------------------------------------
 
    // Condition: a start point is inside the module and outside its dead zone (in strip layers)
    if (!IsPointInsideModule(x, y)) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
 
    // Particle: pdg_code = 0 -> from the MC track of the point refID (as in AddRealPointFull) --
    Int_t pdg = pdg_code;
    if (pdg == 0) {
        FairRunSim* run = FairRunSim::Instance();
        BmnGemStripDigitizer* digi = run ? (BmnGemStripDigitizer*)run->GetTask("GemDigitizer") : nullptr;
        TClonesArray* points = digi ? digi->GetMCPoints() : nullptr;
        TClonesArray* mctracks = digi ? digi->GetMCTracks() : nullptr;
 
        CbmMCTrack* mcTr = nullptr;
        if (points && mctracks && refID >= 0 && refID < points->GetEntriesFast()) {
            CbmStsPoint* point = (CbmStsPoint*)points->UncheckedAt(refID);
            if (point && point->GetTrackID() >= 0 && point->GetTrackID() < mctracks->GetEntriesFast())
                mcTr = (CbmMCTrack*)mctracks->UncheckedAt(point->GetTrackID());
        }
 
        if (mcTr) {
            pdg = mcTr->GetPdgCode();
        } else {
            static std::atomic<bool> no_mc_reported{false};
            if (!no_mc_reported.exchange(true))
                LOG(warn) << "BmnGemStripModule::" << __func__
                          << "(): pdg_code = 0 and the MC track is not available (no GemDigitizer or MC point), "
                             "the He4 mass is used (this warning is printed once)";
        }
    }
 
    TVector3 mom(px, py, pz);
    TParticlePDG* part = (pdg != 0) ? TDatabasePDG::Instance()->GetParticle(pdg) : nullptr;
 
    Double_t mnoc = 1.0, mass = 3.7283999;   // He4
    if (part)
        mass = part->Mass();
    if (mass > 0.0001) {
        TLorentzVector lorv;
        lorv.SetVectM(mom, mass);
        Double_t beta = lorv.Beta();
        Double_t gamma = lorv.Gamma();
        Double_t charge = (part) ? part->Charge() / 3 : 1;
        mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
        if (mnoc < 0.5)
            return kFALSE;
        if (MaxClustersPerCm > 0.0 && mnoc > MaxClustersPerCm)
            mnoc = MaxClustersPerCm;
    }
    //--------------------------------------------------------------------------
 
    // Distance from entry point (x,y,z) to exit point (along z-axis)
    Double_t z_distance_from_in_to_out;
 
    // Condition: track direction along z-axis
    if (pz > 0.0) {
        z_distance_from_in_to_out = ModuleThickness - (z - ZStartModulePosition);
    } else {
        z_distance_from_in_to_out = z - ZStartModulePosition;
    }
 
    // Condition: track distance along z-axis must be not zero
    if (z_distance_from_in_to_out <= 0.0) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") has a track with an incorrect length";
        return false;
    }
 
    // Scale coefficient of the track vector (the path from entry point to exit point)
    Double_t vector_coeff = fabs(z_distance_from_in_to_out / pz);
 
    // Components of the track vector
    Double_t x_vec_from_in_to_out = vector_coeff * px;
    Double_t y_vec_from_in_to_out = vector_coeff * py;
    Double_t z_vec_from_in_to_out = vector_coeff * pz;
 
    // Exit point coordinates (x_out, y_out, z_out)
    Double_t x_out = x + x_vec_from_in_to_out;
    Double_t y_out = y + y_vec_from_in_to_out;
    Double_t z_out = z + z_vec_from_in_to_out;
 
    // Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
    // Condition: x-coordinate of the exit point is inside the module
    if (x_out < XMinModule || x_out > XMaxModule) {
        if (x_out < XMinModule)
            x_out = XMinModule;
        if (x_out > XMaxModule)
            x_out = XMaxModule;
 
        x_vec_from_in_to_out = x_out - x;
        vector_coeff = x_vec_from_in_to_out / px;
 
        y_vec_from_in_to_out = vector_coeff * py;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        y_out = y + y_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    // Condition: y-coordinate of the exit point is inside the module
    if (y_out < YMinModule || y_out > YMaxModule) {
        if (y_out < YMinModule)
            y_out = YMinModule;
        if (y_out > YMaxModule)
            y_out = YMaxModule;
 
        y_vec_from_in_to_out = y_out - y;
        vector_coeff = y_vec_from_in_to_out / py;
 
        x_vec_from_in_to_out = vector_coeff * px;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        x_out = x + x_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    //--------------------------------------------------------------------------
 
    // Center point of the ANN strip window: the track at the middle of the drift gap (computed always, also with
    // onlyStandard: it is stored with the point, see GetRealPointXcenter(), for the training samples)
    Double_t x_c = x, y_c = y;
    GetANNWindowCenterPoint(x, y, z, x_out, y_out, z_out, x_c, y_c);
 
    // Common track length
    Double_t track_length =
        std::sqrt((x_vec_from_in_to_out) * (x_vec_from_in_to_out) + (y_vec_from_in_to_out) * (y_vec_from_in_to_out)
                  + (z_vec_from_in_to_out) * (z_vec_from_in_to_out));
 
    // Which method is used for this point: standard (electrons -> strips) or ANN (model -> strips)
    // (without a suitable model or with onlyStandard the standard method is used)
    Double_t r_vec = TMath::Sqrt(px * px + py * py + pz * pz);
    Double_t theta_deg_abs = TMath::ACos(TMath::Abs(pz) / r_vec) * TMath::RadToDeg();
    const Bool_t useANN = (theta_deg_abs < maxThetaAngle) && mlp_model;
 
    // ANN with external segments: the avalanche is not needed
    const Bool_t use_external_segs = useANN && !segs.empty();
    if (use_external_segs && segs.size() < static_cast<size_t>(n_segs_features)) {
        LOG(error) << "BmnGemStripModule::" << __func__ << "(): segment vector size = " << segs.size()
                   << ", the model requires " << n_segs_features;
        return false;
    }
 
    // Each level - distance to the readout plane
    Double_t level1 = InductionGapThickness;
    Double_t level2 = InductionGapThickness + SecondTransferGapThickness;
    Double_t level3 = InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness;
 
    // Mean electron shift along x-axis (under the influence of the Lorentz force)
    Double_t xmean;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
 
    // Sigma electron smearing
    Double_t sigma;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3+ p4*x^4 + p5*x^5)
 
    // Reading the medium configuration -------------------------------------
    Double_t p0_xmean = BmnGemStripMedium::GetInstance().p0_xmean;
    Double_t p1_xmean = BmnGemStripMedium::GetInstance().p1_xmean;
    Double_t p2_xmean = BmnGemStripMedium::GetInstance().p2_xmean;
    Double_t p3_xmean = BmnGemStripMedium::GetInstance().p3_xmean;
 
    Double_t p0_sigma = BmnGemStripMedium::GetInstance().p0_sigma;
    Double_t p1_sigma = BmnGemStripMedium::GetInstance().p1_sigma;
    Double_t p2_sigma = BmnGemStripMedium::GetInstance().p2_sigma;
    Double_t p3_sigma = BmnGemStripMedium::GetInstance().p3_sigma;
    Double_t p4_sigma = BmnGemStripMedium::GetInstance().p4_sigma;
    Double_t p5_sigma = BmnGemStripMedium::GetInstance().p5_sigma;
    //----------------------------------------------------------------------
 
    // Initialization of a strip cluster in each strip layer
    Int_t n_strip_layers = StripLayers.size();   // number of layers in the module
    vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
    // Segmented charge distribution of this point (always computed when the avalanche is simulated)
    vector<Int_t> seg_local(n_track_segments, 0);
 
    // Standard method: per-layer strip accumulators -------------------------
    // Signals are accumulated per strip and added to the cluster with one AddStrip() per strip:
    // AddStrip() sums the signals and keeps the strips in ascending order, the signals are integer
    // numbers of electrons -> the cluster is bitwise the same as with one AddStrip() per electron.
    vector<vector<Double_t>> strip_acc;
    vector<Int_t> acc_min, acc_max;   // range of the touched strips in each layer
    vector<char> layer_vertical;      // layer invariant: vertical strips (|angle| < 1 deg)
    if (!useANN) {
        strip_acc.resize(n_strip_layers);
        acc_min.assign(n_strip_layers, std::numeric_limits<Int_t>::max());
        acc_max.assign(n_strip_layers, -1);
        layer_vertical.assign(n_strip_layers, 0);
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            strip_acc[ilayer].assign(std::max(0, StripLayers[ilayer].GetNStrips()), 0.0);
            layer_vertical[ilayer] = TMath::Abs(StripLayers[ilayer].GetAngleDeg()) < 1.0;
        }
    }
    //--------------------------------------------------------------------------
 
    if (!use_external_segs) {
 
        Double_t current_length = 0.0;         // traversed length (counter)
        Double_t current_length_ratio = 0.0;   // ratio of the traversed length to common track length (counter)
        Double_t current_step = 0.0;           // current distance between two collision points
 
        // Collection of collision points
        std::vector<CollPoint> collision_points;
 
        // Reading the MCD value (mean collision distance = mean free flight path) [cm]
        MCD = BmnGemStripMedium::GetInstance().MCD;
        MCD = 1 / mnoc;   // AZ-110124
 
        while (current_length < track_length) {
 
            current_step = gRandom->Exp(MCD);
            current_length += current_step;
 
            if (current_length > track_length)
                break;
 
            current_length_ratio = current_length / track_length;
 
            Double_t current_x = x + current_length_ratio * x_vec_from_in_to_out;
            Double_t current_y = y + current_length_ratio * y_vec_from_in_to_out;
            Double_t current_z = z + current_length_ratio * z_vec_from_in_to_out;
 
            collision_points.push_back(CollPoint(current_x, current_y, current_z));
        }
 
        // Electron avalanche producing for each collision point at the track
        // (in the ANN mode only the segmented charge distribution is filled here)
        for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
 
            Double_t xcoll = collision_points[icoll].x;
            Double_t ycoll = collision_points[icoll].y;
            Double_t zcoll = collision_points[icoll].z;
 
            Double_t zdist;   // current z-distance to the readout
 
            // Find z-distance to the readout depending on the electron drift direction
            if (ElectronDriftDirection == ForwardZAxisEDrift) {
                zdist = (ModuleThickness + ZStartModulePosition) - zcoll;
            } else {
                zdist = zcoll - ZStartModulePosition;
            }
            //------------------------------------------------------------------
 
            // polynomials
            xmean = p0_xmean + p1_xmean * zdist + p2_xmean * zdist * zdist + p3_xmean * zdist * zdist * zdist;
            sigma = p0_sigma + p1_sigma * zdist + p2_sigma * zdist * zdist + p3_sigma * zdist * zdist * zdist
                    + p4_sigma * zdist * zdist * zdist * zdist + p5_sigma * zdist * zdist * zdist * zdist * zdist;
 
            // Number of electrons in the current collision
            Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
            if (n_electrons_cluster < 1)
                n_electrons_cluster = 1;   // min
            if (n_electrons_cluster > 6)
                n_electrons_cluster = 6;   // max
 
            // Filling charge distribution -------------------------------------
            Double_t dist_from_start = 0.0;
            if (pz > 0.0) {
                dist_from_start =
                    std::sqrt((xcoll - x) * (xcoll - x) + (ycoll - y) * (ycoll - y) + (zcoll - z) * (zcoll - z));
            } else {
                dist_from_start = std::sqrt((xcoll - x_out) * (xcoll - x_out) + (ycoll - y_out) * (ycoll - y_out)
                                            + (zcoll - z_out) * (zcoll - z_out));
            }
            Int_t segment = (Int_t)(dist_from_start / track_length * n_track_segments);
            if (segment >= n_track_segments)
                segment = n_track_segments - 1;
            //------------------------------------------------------------------
 
            for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
 
                // Electron gain in each GEM cascade (three draws always, as in the other methods)
                Double_t gain_gem1 = gRandom->Exp(15);
                Double_t gain_gem2 = gRandom->Exp(15);
                Double_t gain_gem3 = gRandom->Exp(15);
 
                if (gain_gem1 < 1.0)
                    gain_gem1 = 1.0;
                if (gain_gem2 < 1.0)
                    gain_gem2 = 1.0;
                if (gain_gem3 < 1.0)
                    gain_gem3 = 1.0;
 
                int total_gain = 0;
 
                if (zdist < level1) {
                    total_gain = 1.0;
                }
                if (zdist >= level1 && zdist < level2) {
                    total_gain = gain_gem3;
                }
                if (zdist >= level2 && zdist < level3) {
                    total_gain = gain_gem3 * gain_gem2;
                }
                if (zdist >= level3) {
                    total_gain = gain_gem3 * gain_gem2 * gain_gem1;
                }
 
                seg_local[segment] += total_gain;
 
                // ANN mode: electrons are not projected on the readout, strips are taken from the model
                if (useANN)
                    continue;
 
                // Standard method ---------------------------------------------
                // Projection of the current electron on the readout (x,y-coordinates)
                double x_readout, y_readout;
 
                for (int igain = 0; igain < total_gain; ++igain) {
 
                    // x-shift of the electon depends on the electron drift direction
                    if (ElectronDriftDirection == ForwardZAxisEDrift) {
                        x_readout = gRandom->Gaus(xcoll + xmean, sigma);
                    } else {
                        x_readout = gRandom->Gaus(xcoll - xmean, sigma);
                    }
 
                    y_readout = gRandom->Gaus(ycoll, sigma);
 
                    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                        // Condition: end electron position must be inside the layer
                        if (!StripLayers[ilayer].IsPointInsideStripLayer(x_readout, y_readout))
                            continue;
 
                        Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_readout, y_readout);
                        // AZ-100124 - account for strip width
                        Double_t dstrip = strip_pos - (int)strip_pos;
                        if (layer_vertical[ilayer]) {
                            // Vertical strips
                            if (dstrip >= 0.4 && dstrip < 0.6)
                                continue;
                        } else {
                            // Rotated strips
                            if (dstrip >= 0.3 && dstrip < 0.7)
                                continue;
                        }
                        // AZ
 
                        const Int_t strip_num = (Int_t)strip_pos;
                        if (strip_num >= 0 && strip_num < (Int_t)strip_acc[ilayer].size()) {
                            strip_acc[ilayer][strip_num] += 1.0;
                            if (strip_num < acc_min[ilayer])
                                acc_min[ilayer] = strip_num;
                            if (strip_num > acc_max[ilayer])
                                acc_max[ilayer] = strip_num;
                        } else {
                            // out of the accumulator range: directly, as in the other methods
                            cluster_layers[ilayer].AddStrip(strip_num, 1.0);
                        }
                    }
 
                    // Electron points on the readout (generation data) ---------
                    if (FillGenerationData) {
                        XElectronPos.push_back(x_readout);
                        YElectronPos.push_back(y_readout);
                        ElectronSignal.push_back(1.0);
                    }
                    //----------------------------------------------------------
                }
            }
        }
    }
 
    // Standard method: accumulated signals -> clusters (one AddStrip per strip)
    if (!useANN) {
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            for (Int_t istrip = acc_min[ilayer]; istrip <= acc_max[ilayer]; ++istrip) {
                if (strip_acc[ilayer][istrip] > 0.0)
                    cluster_layers[ilayer].AddStrip(istrip, strip_acc[ilayer][istrip]);
            }
        }
    }
    //----------------------------------------------------------------------
 
    // ANN generation (MLP, own inference) ---------------------------------
    if (useANN) {
 
        // Segments for the model: external or from the avalanche of this point
        const vector<Int_t>& seg_used = use_external_segs ? segs : seg_local;
 
        // Normalization (px,py,pz), where r_vec must be equal to 1
        Double_t px_norm = px / r_vec;
        Double_t py_norm = py / r_vec;
        Double_t pz_norm = pz / r_vec;
 
        // Inversion for z-reverse direction particles: the same track segment passed from (x_out, y_out)
        if (pz < 0.0) {
            px_norm = -px_norm;
            py_norm = -py_norm;
            pz_norm = -pz_norm;
        }
 
        Double_t drift_direction = 0;   // 0 is ForwardZAxisEDrift
        if (ElectronDriftDirection == BackwardZAxisEDrift)
            drift_direction = 1;
 
        // Magnetic field: given explicitly or from the current medium configuration
        const Double_t bfield_used = (bfield < 0.0) ? BmnGemStripMedium::GetInstance().BField : bfield;
 
        // Normalize seg. vector: [0, max] -> [0, 1] (max over the whole vector)
        Double_t max_segment_val = *(std::max_element(seg_used.begin(), seg_used.end()));
 
        // Total charge in segment vector
        Double_t total_charge = 0.0;
        for (size_t i = 0; i < seg_used.size(); ++i)
            total_charge += seg_used[i];
 
        // Window center point (x_c, y_c): computed above with GetANNWindowCenterPoint()
 
        // Common part of the model input: [px, py, pz, drift, bfield, s0..s{n_segs_features-1}]
        float common_input[n_common_inputs];
        common_input[0] = static_cast<float>(px_norm);
        common_input[1] = static_cast<float>(py_norm);
        common_input[2] = static_cast<float>(pz_norm);
        common_input[3] = static_cast<float>(drift_direction);
        common_input[4] = static_cast<float>(bfield_used);
        for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
            common_input[5 + iseg] =
                (max_segment_val > 0.0) ? static_cast<float>(seg_used[iseg] / max_segment_val) : 0.0f;
 
        if (total_charge > 0.0) {
            // Strip model: a window of n_ann_strips strips around the window center strip in each layer
            const Int_t center_strip_index = n_ann_strips / 2;   // center strip of the window
 
            // Layers with the MC point inside: their part of the input [strip_layer, strip_order, center_frac]
            std::vector<Int_t> ann_layers;          // layer index
            std::vector<Int_t> ann_center_strips;   // MC strip (window center)
            std::vector<Int_t> ann_n_strips;        // number of strips in the layer
            std::vector<float> ann_layer_inputs;    // [n_ann_layers][n_layer_inputs]
            ann_layers.reserve(n_strip_layers);
            ann_center_strips.reserve(n_strip_layers);
            ann_n_strips.reserve(n_strip_layers);
            ann_layer_inputs.reserve(n_strip_layers * n_layer_inputs);
 
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
                // The window is built around the center point (middle of the drift gap): it must be inside this layer
                if (!StripLayers[ilayer].IsPointInsideStripLayer(x_c, y_c))
                    continue;
 
                Int_t n_strips_in_layer = StripLayers[ilayer].GetNStrips();
                if (n_strips_in_layer <= 0)
                    continue;
 
                // Center point position in strips: integer part - window center strip, fractional part - center_frac
                Double_t center_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_c, y_c);
                Int_t center_strip = (Int_t)TMath::Floor(center_pos);
                Double_t center_frac = center_pos - center_strip;
 
                Double_t strip_layer = static_cast<Double_t>(StripLayers[ilayer].GetType());   // 0 - lower, 1 - upper
                Double_t strip_order =
                    static_cast<Double_t>(StripLayers[ilayer].GetStripNumberingOrder());   // 0 - LeftToRight
 
                ann_layers.push_back(ilayer);
                ann_center_strips.push_back(center_strip);
                ann_n_strips.push_back(n_strips_in_layer);
                ann_layer_inputs.push_back(static_cast<float>(strip_layer));
                ann_layer_inputs.push_back(static_cast<float>(strip_order));
                ann_layer_inputs.push_back(static_cast<float>(center_frac));
            }
 
            // Strip signals for all these layers in one call: n_ann_strips values in [0, 1] (max ~ 1) per layer
            Int_t n_ann_layers = ann_layers.size();
            std::vector<float> ann_strip_values(static_cast<size_t>(n_ann_layers) * n_ann_strips);
            if (n_ann_layers > 0
                && !mlp_model->PredictSplit(common_input, n_common_inputs, ann_layer_inputs.data(), n_ann_layers,
                                            ann_strip_values.data()))
            {
                LOG(error) << "BmnGemStripModule::" << __func__ << "(): MLP prediction failed";
                n_ann_layers = 0;
            }
 
            for (Int_t iann = 0; iann < n_ann_layers; ++iann) {
                const Int_t ilayer = ann_layers[iann];
                const Int_t center_strip = ann_center_strips[iann];
                const Int_t n_strips_in_layer = ann_n_strips[iann];
                const float* strip_values = ann_strip_values.data() + static_cast<size_t>(iann) * n_ann_strips;
 
                // Sum of the model output over all strips of the window
                Double_t sum_norm = 0.0;
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip)
                    sum_norm += strip_values[istrip];
 
                // Scale coeff for total signal distribution (separate for lower and upper layers)
                const Double_t ann_correction_lower = 0.8;   // Set to 1.0 to disable correction
                const Double_t ann_correction_upper = 0.6;   // Set to 1.0 to disable correction
                const Double_t ann_correction =
                    (StripLayers[ilayer].GetType() == LowerStripLayer) ? ann_correction_lower : ann_correction_upper;
 
                Double_t signal_scale = 0.0;
                if (sum_norm > 0.0)
                    signal_scale = (total_charge / sum_norm) * ann_correction;
 
                Double_t threshold = 0.0001 * total_charge;
 
                // Strips are written directly: window index -> strip number
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip) {
                    Double_t content = strip_values[istrip] * signal_scale;
                    if (content <= threshold)
                        continue;
 
                    Int_t strip_num = center_strip - center_strip_index + istrip;
                    if (strip_num < 0 || strip_num >= n_strips_in_layer)
                        continue;   // the window goes out of the layer
 
                    cluster_layers[ilayer].AddStrip(strip_num, content);
                }
            }
        }
    }
    //----------------------------------------------------------------------
 
    // Segmented charge distribution (generation data) ------------------------
    // (seg_used above is not used after this point)
    if (FillGenerationData) {
        if (use_external_segs)
            SegmentedChargeDistribution = std::move(segs);
        else
            SegmentedChargeDistribution = std::move(seg_local);
    }
    //----------------------------------------------------------------------
 
    // Calculate cluster parameters --------------------------------------------
    vector<Double_t> cluster_mean_position(n_strip_layers, 0.0);
    vector<Double_t> cluster_total_signal(n_strip_layers, 0.0);
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
        if (cluster_layers[ilayer].GetClusterSize() > 0) {
 
            for (int i = 0; i < cluster_layers[ilayer].GetClusterSize(); ++i) {
                cluster_mean_position[ilayer] += (cluster_layers[ilayer].Strips[i] + 0.5)
                                                 * cluster_layers[ilayer].Signals[i];   // as sum of all positions
                cluster_total_signal[ilayer] += cluster_layers[ilayer].Signals[i];
            }
            cluster_mean_position[ilayer] /= cluster_total_signal[ilayer];
        }
 
        cluster_layers[ilayer].MeanPosition = cluster_mean_position[ilayer];   // mean cluster position
        cluster_layers[ilayer].TotalSignal = cluster_total_signal[ilayer];     // total signal of the cluster
    }
 
    if (ElectronDriftDirection == ForwardZAxisEDrift) {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        }
    } else {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        cluster_layers[ilayer].PositionResidual =
            cluster_layers[ilayer].MeanPosition - cluster_layers[ilayer].OriginPosition;
        cluster_layers[ilayer].SetType(StripLayers[ilayer].GetType());
    }
    //----------------------------------------------------------------------
 
    // Added clusters (generation data) ----------------------------------------
    if (FillGenerationData) {
        if (AddedClusters.size() == 0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                AddedClusters.push_back(vector<StripCluster>());
            }
        }
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
        }
    }
    //------------------------------------------------------------------------------
 
    // Add the correct clusters to the strip layers -------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
        }
    }
    //----------------------------------------------------------------------
 
    // Fill strip matches ---------------------------------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                    refID);
        }
    }
    //----------------------------------------------------------------------
 
    RealPointsX.push_back(x);
    RealPointsY.push_back(y);
    RealPointsXout.push_back(x_out);
    RealPointsYout.push_back(y_out);
    RealPointsMC.push_back(refID);
    RealPointsXcenter.push_back(x_c);   // ANN window center point
    RealPointsYcenter.push_back(y_c);
 
    return true;
}
//------------------------------------------------------------------------------
 
// The same as AddRealPointFull_ANNStripGenerator_MLP, the standard method is parallel (std::thread);
Bool_t BmnGemStripModule::AddRealPointFull_ANNStripGenerator_MLP_Parallel(Double_t x,
                                                                          Double_t y,
                                                                          Double_t z,
                                                                          Double_t px,
                                                                          Double_t py,
                                                                          Double_t pz,
                                                                          Double_t signal,
                                                                          Int_t refID,
                                                                          Int_t pdg_code,
                                                                          Double_t bfield,
                                                                          Bool_t onlyStandard,
                                                                          vector<Int_t> segs)
{
    const Double_t maxThetaAngle = kANNMaxThetaAngle;   // max |theta| (deg) of the track for the ANN
 
    // Model format: [px, py, pz, drift, bfield, s0..s63, strip_layer, strip_order, center_frac] -> window of N
    // strips (N odd, from the model) centered at the middle of the drift gap (GetANNWindowCenterPoint)
    const Int_t n_track_segments = 64;                   // track segmentation (segmented charge distribution)
    const Int_t n_segs_features = n_track_segments;      // segments passed to the model
    const Int_t n_common_inputs = 5 + n_segs_features;   // the same for all strip layers of the point
    const Int_t n_layer_inputs = 3;                      // [strip_layer, strip_order, center_frac]
    Int_t n_ann_strips = 0;   // model output: strip window size (odd), taken from the model below
 
    // Test outputs are cleared always (cheap if they are empty)
    ResetElectronPointsAndClusters();
 
    // MLP model: shared between modules, checked once for each model ---------
    std::shared_ptr<const BmnGemStripMLPModel> mlp_holder;
    const BmnGemStripMLPModel* mlp_model = nullptr;
    if (!onlyStandard) {
        mlp_holder = BmnGemStripMLPModel::GetShared();
        mlp_model = mlp_holder.get();
 
        static std::mutex checked_mutex;
        static std::map<const BmnGemStripMLPModel*, std::pair<std::shared_ptr<const BmnGemStripMLPModel>, std::string>>
            checked;   // model -> (holder, why it is not suitable: empty if suitable)
 
        std::lock_guard<std::mutex> lock(checked_mutex);
        if (!mlp_model) {
            LOG(warn) << "BmnGemStripModule::" << __func__
                      << "(): MLP model is not available, the standard method is used";
        } else {
            auto it = checked.find(mlp_model);
            if (it == checked.end()) {
                std::vector<std::string> expected = {"px", "py", "pz", "drift", "bfield"};
                for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
                    expected.push_back("s" + std::to_string(iseg));
                expected.insert(expected.end(), {"strip_layer", "strip_order", "center_frac"});
 
                std::string why;
                Bool_t ok = mlp_model->CheckFeatureNames(expected, &why);
                if (ok) {
                    why = CheckANNStripWindowSize(mlp_model->GetOutputDim());
                    ok = why.empty();
                }
                if (ok) {
                    why.clear();
                    LOG(info) << "BmnGemStripModule::" << __func__ << "(): MLP model " << mlp_model->GetFilePath()
                              << " is used (ANN for |theta| < " << maxThetaAngle << " deg, window "
                              << mlp_model->GetOutputDim() << " strips)";
                } else if (why.empty())
                    why = "unknown reason";
                it = checked.emplace(mlp_model, std::make_pair(mlp_holder, why)).first;
            }
            if (!it->second.second.empty()) {
                LOG(warn) << "BmnGemStripModule::" << __func__ << "(): the model " << mlp_model->GetFilePath()
                          << " is not suitable (" << it->second.second << "), the standard method is used";
                mlp_model = nullptr;
            }
        }
    }
    if (mlp_model)
        n_ann_strips = mlp_model->GetOutputDim();   // window size of the model (checked above: odd, >= 3)
    //--------------------------------------------------------------------------
 
    // Condition: a start point is inside the module and outside its dead zone (in strip layers)
    if (!IsPointInsideModule(x, y)) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
 
    // Particle: pdg_code = 0 -> from the MC track of the point refID (as in AddRealPointFull) --
    Int_t pdg = pdg_code;
    if (pdg == 0) {
        FairRunSim* run = FairRunSim::Instance();
        BmnGemStripDigitizer* digi = run ? (BmnGemStripDigitizer*)run->GetTask("GemDigitizer") : nullptr;
        TClonesArray* points = digi ? digi->GetMCPoints() : nullptr;
        TClonesArray* mctracks = digi ? digi->GetMCTracks() : nullptr;
 
        CbmMCTrack* mcTr = nullptr;
        if (points && mctracks && refID >= 0 && refID < points->GetEntriesFast()) {
            CbmStsPoint* point = (CbmStsPoint*)points->UncheckedAt(refID);
            if (point && point->GetTrackID() >= 0 && point->GetTrackID() < mctracks->GetEntriesFast())
                mcTr = (CbmMCTrack*)mctracks->UncheckedAt(point->GetTrackID());
        }
 
        if (mcTr) {
            pdg = mcTr->GetPdgCode();
        } else {
            static std::atomic<bool> no_mc_reported{false};
            if (!no_mc_reported.exchange(true))
                LOG(warn) << "BmnGemStripModule::" << __func__
                          << "(): pdg_code = 0 and the MC track is not available (no GemDigitizer or MC point), "
                             "the He4 mass is used (this warning is printed once)";
        }
    }
 
    TVector3 mom(px, py, pz);
    TParticlePDG* part = (pdg != 0) ? TDatabasePDG::Instance()->GetParticle(pdg) : nullptr;
 
    Double_t mnoc = 1.0, mass = 3.7283999;   // He4
    if (part)
        mass = part->Mass();
    if (mass > 0.0001) {
        TLorentzVector lorv;
        lorv.SetVectM(mom, mass);
        Double_t beta = lorv.Beta();
        Double_t gamma = lorv.Gamma();
        Double_t charge = (part) ? part->Charge() / 3 : 1;
        mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
        if (mnoc < 0.5)
            return kFALSE;
        if (MaxClustersPerCm > 0.0 && mnoc > MaxClustersPerCm)
            mnoc = MaxClustersPerCm;
    }
    //--------------------------------------------------------------------------
 
    // Distance from entry point (x,y,z) to exit point (along z-axis)
    Double_t z_distance_from_in_to_out;
 
    // Condition: track direction along z-axis
    if (pz > 0.0) {
        z_distance_from_in_to_out = ModuleThickness - (z - ZStartModulePosition);
    } else {
        z_distance_from_in_to_out = z - ZStartModulePosition;
    }
 
    // Condition: track distance along z-axis must be not zero
    if (z_distance_from_in_to_out <= 0.0) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") has a track with an incorrect length";
        return false;
    }
 
    // Scale coefficient of the track vector (the path from entry point to exit point)
    Double_t vector_coeff = fabs(z_distance_from_in_to_out / pz);
 
    // Components of the track vector
    Double_t x_vec_from_in_to_out = vector_coeff * px;
    Double_t y_vec_from_in_to_out = vector_coeff * py;
    Double_t z_vec_from_in_to_out = vector_coeff * pz;
 
    // Exit point coordinates (x_out, y_out, z_out)
    Double_t x_out = x + x_vec_from_in_to_out;
    Double_t y_out = y + y_vec_from_in_to_out;
    Double_t z_out = z + z_vec_from_in_to_out;
 
    // Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
    // Condition: x-coordinate of the exit point is inside the module
    if (x_out < XMinModule || x_out > XMaxModule) {
        if (x_out < XMinModule)
            x_out = XMinModule;
        if (x_out > XMaxModule)
            x_out = XMaxModule;
 
        x_vec_from_in_to_out = x_out - x;
        vector_coeff = x_vec_from_in_to_out / px;
 
        y_vec_from_in_to_out = vector_coeff * py;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        y_out = y + y_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    // Condition: y-coordinate of the exit point is inside the module
    if (y_out < YMinModule || y_out > YMaxModule) {
        if (y_out < YMinModule)
            y_out = YMinModule;
        if (y_out > YMaxModule)
            y_out = YMaxModule;
 
        y_vec_from_in_to_out = y_out - y;
        vector_coeff = y_vec_from_in_to_out / py;
 
        x_vec_from_in_to_out = vector_coeff * px;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        x_out = x + x_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    //--------------------------------------------------------------------------
 
    // Center point of the ANN strip window: the track at the middle of the drift gap (computed always, also with
    // onlyStandard: it is stored with the point, see GetRealPointXcenter(), for the training samples)
    Double_t x_c = x, y_c = y;
    GetANNWindowCenterPoint(x, y, z, x_out, y_out, z_out, x_c, y_c);
 
    // Common track length
    Double_t track_length =
        std::sqrt((x_vec_from_in_to_out) * (x_vec_from_in_to_out) + (y_vec_from_in_to_out) * (y_vec_from_in_to_out)
                  + (z_vec_from_in_to_out) * (z_vec_from_in_to_out));
 
    // Which method is used for this point: standard (electrons -> strips) or ANN (model -> strips)
    // (without a suitable model or with onlyStandard the standard method is used)
    Double_t r_vec = TMath::Sqrt(px * px + py * py + pz * pz);
    Double_t theta_deg_abs = TMath::ACos(TMath::Abs(pz) / r_vec) * TMath::RadToDeg();
    const Bool_t useANN = (theta_deg_abs < maxThetaAngle) && mlp_model;
 
    // ANN with external segments: the avalanche is not needed
    const Bool_t use_external_segs = useANN && !segs.empty();
    if (use_external_segs && segs.size() < static_cast<size_t>(n_segs_features)) {
        LOG(error) << "BmnGemStripModule::" << __func__ << "(): segment vector size = " << segs.size()
                   << ", the model requires " << n_segs_features;
        return false;
    }
 
    // Each level - distance to the readout plane
    Double_t level1 = InductionGapThickness;
    Double_t level2 = InductionGapThickness + SecondTransferGapThickness;
    Double_t level3 = InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness;
 
    // Mean electron shift along x-axis (under the influence of the Lorentz force)
    Double_t xmean;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
 
    // Sigma electron smearing
    Double_t sigma;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3+ p4*x^4 + p5*x^5)
 
    // Reading the medium configuration -------------------------------------
    Double_t p0_xmean = BmnGemStripMedium::GetInstance().p0_xmean;
    Double_t p1_xmean = BmnGemStripMedium::GetInstance().p1_xmean;
    Double_t p2_xmean = BmnGemStripMedium::GetInstance().p2_xmean;
    Double_t p3_xmean = BmnGemStripMedium::GetInstance().p3_xmean;
 
    Double_t p0_sigma = BmnGemStripMedium::GetInstance().p0_sigma;
    Double_t p1_sigma = BmnGemStripMedium::GetInstance().p1_sigma;
    Double_t p2_sigma = BmnGemStripMedium::GetInstance().p2_sigma;
    Double_t p3_sigma = BmnGemStripMedium::GetInstance().p3_sigma;
    Double_t p4_sigma = BmnGemStripMedium::GetInstance().p4_sigma;
    Double_t p5_sigma = BmnGemStripMedium::GetInstance().p5_sigma;
    //----------------------------------------------------------------------
 
    // Initialization of a strip cluster in each strip layer
    Int_t n_strip_layers = StripLayers.size();   // number of layers in the module
    vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
    // Segmented charge distribution of this point (always computed when the avalanche is simulated)
    vector<Int_t> seg_local(n_track_segments, 0);
 
    // [PARALLEL-STD] Standard method: parallel processing of the avalanche electrons -----
    // true  - electron (x,y) are taken from gRandom sequentially, as in AddRealPointFull_ANNStripGenerator_MLP:
    //         bitwise identical result, only the strip assignment is parallel;
    // false - electron (x,y) are generated in parallel by per-task generators (seeded from gRandom):
    //         faster, statistically equivalent to the sequential methods, bitwise identical to
    //         AddRealPointFull_ANNStripGenerator_Fdeep_Parallel (the same scheme of the random numbers).
    constexpr Bool_t kExactGRandomSequence = false;
    constexpr Long64_t kElectronsPerTask = 4096;           // electrons per parallel task
    constexpr Long64_t kMinElectronsForParallel = 20000;   // below this - sequential execution
 
    struct AvalancheTask   // used only if !kExactGRandomSequence
    {
        Double_t xcenter;   // mean x on the readout (with the Lorentz shift)
        Double_t ycenter;   // mean y on the readout
        Double_t sigma;     // smearing
        Long64_t offset;    // index of the first electron of the task (from first_electron)
        Long64_t n;         // number of electrons in the task
    };
    std::vector<AvalancheTask> aval_tasks;
    Long64_t n_aval_electrons = 0;
    std::vector<Double_t> exact_x, exact_y;   // electron positions, used only if kExactGRandomSequence
 
    // Persistent thread pool: the threads are created once (on the first call) and wait for work,
    // so no threads are created/joined in each call. Lives until the end of the program.
    // Work items are taken dynamically (atomic counter) -> load balancing for avalanches of different size.
    struct ThreadPool
    {
        std::vector<std::thread> workers;
        std::mutex mtx;                     // protects the job description and the counters below
        std::mutex run_mtx;                 // one job at a time (in case of calls from several threads)
        std::condition_variable cv_start;   // workers: a new job is available
        std::condition_variable cv_done;    // caller: all workers have finished the job
        const std::function<void(size_t)>* job = nullptr;
        size_t n_items = 0;
        std::atomic<size_t> next{0};   // next work item
        size_t generation = 0;         // job counter
        size_t n_busy = 0;             // workers that have not finished the current job yet
        bool stop = false;
 
        explicit ThreadPool(size_t n_workers)
        {
            workers.reserve(n_workers);
            for (size_t i = 0; i < n_workers; ++i)
                workers.emplace_back([this] { Loop(); });
        }
 
        ~ThreadPool()
        {
            {
                std::lock_guard<std::mutex> lk(mtx);
                stop = true;
            }
            cv_start.notify_all();
            for (auto& w : workers)
                w.join();
        }
 
        void Work()
        {
            for (size_t i = next.fetch_add(1); i < n_items; i = next.fetch_add(1))
                (*job)(i);
        }
 
        void Loop()
        {
            size_t seen = 0;
            for (;;) {
                {
                    std::unique_lock<std::mutex> lk(mtx);
                    cv_start.wait(lk, [&] { return stop || generation != seen; });
                    if (stop)
                        return;
                    seen = generation;
                }
                Work();
                {
                    std::lock_guard<std::mutex> lk(mtx);
                    if (--n_busy == 0)
                        cv_done.notify_one();
                }
            }
        }
 
        // Runs func(i) for i in [0, n) on the pool threads and the calling thread
        void Run(size_t n, const std::function<void(size_t)>& func)
        {
            if (n == 0)
                return;
            if (workers.empty() || n == 1) {
                for (size_t i = 0; i < n; ++i)
                    func(i);
                return;
            }
 
            std::lock_guard<std::mutex> run_lock(run_mtx);
            {
                std::lock_guard<std::mutex> lk(mtx);
                job = &func;
                n_items = n;
                next.store(0);
                n_busy = workers.size();
                ++generation;
            }
            cv_start.notify_all();
 
            Work();   // the calling thread works too
 
            std::unique_lock<std::mutex> lk(mtx);
            cv_done.wait(lk, [&] { return n_busy == 0; });
            job = nullptr;
        }
    };
 
    // hardware threads - 1 workers (the calling thread is the last one)
    static ThreadPool pool(std::max(1u, std::thread::hardware_concurrency()) - 1);
 
    auto parallel_for = [](size_t n, const std::function<void(size_t)>& func) { pool.Run(n, func); };
    //--------------------------------------------------------------------------
 
    if (!use_external_segs) {
 
        Double_t current_length = 0.0;         // traversed length (counter)
        Double_t current_length_ratio = 0.0;   // ratio of the traversed length to common track length (counter)
        Double_t current_step = 0.0;           // current distance between two collision points
 
        // Collection of collision points
        std::vector<CollPoint> collision_points;
 
        // Reading the MCD value (mean collision distance = mean free flight path) [cm]
        MCD = BmnGemStripMedium::GetInstance().MCD;
        MCD = 1 / mnoc;   // AZ-110124
 
        while (current_length < track_length) {
 
            current_step = gRandom->Exp(MCD);
            current_length += current_step;
 
            if (current_length > track_length)
                break;
 
            current_length_ratio = current_length / track_length;
 
            Double_t current_x = x + current_length_ratio * x_vec_from_in_to_out;
            Double_t current_y = y + current_length_ratio * y_vec_from_in_to_out;
            Double_t current_z = z + current_length_ratio * z_vec_from_in_to_out;
 
            collision_points.push_back(CollPoint(current_x, current_y, current_z));
        }
 
        // Electron avalanche producing for each collision point at the track
        // (in the ANN mode only the segmented charge distribution is filled here)
        for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
 
            Double_t xcoll = collision_points[icoll].x;
            Double_t ycoll = collision_points[icoll].y;
            Double_t zcoll = collision_points[icoll].z;
 
            Double_t zdist;   // current z-distance to the readout
 
            // Find z-distance to the readout depending on the electron drift direction
            if (ElectronDriftDirection == ForwardZAxisEDrift) {
                zdist = (ModuleThickness + ZStartModulePosition) - zcoll;
            } else {
                zdist = zcoll - ZStartModulePosition;
            }
            //------------------------------------------------------------------
 
            // polynomials
            xmean = p0_xmean + p1_xmean * zdist + p2_xmean * zdist * zdist + p3_xmean * zdist * zdist * zdist;
            sigma = p0_sigma + p1_sigma * zdist + p2_sigma * zdist * zdist + p3_sigma * zdist * zdist * zdist
                    + p4_sigma * zdist * zdist * zdist * zdist + p5_sigma * zdist * zdist * zdist * zdist * zdist;
 
            // Number of electrons in the current collision
            Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
            if (n_electrons_cluster < 1)
                n_electrons_cluster = 1;   // min
            if (n_electrons_cluster > 6)
                n_electrons_cluster = 6;   // max
 
            // Filling charge distribution -------------------------------------
            Double_t dist_from_start = 0.0;
            if (pz > 0.0) {
                dist_from_start =
                    std::sqrt((xcoll - x) * (xcoll - x) + (ycoll - y) * (ycoll - y) + (zcoll - z) * (zcoll - z));
            } else {
                dist_from_start = std::sqrt((xcoll - x_out) * (xcoll - x_out) + (ycoll - y_out) * (ycoll - y_out)
                                            + (zcoll - z_out) * (zcoll - z_out));
            }
            Int_t segment = (Int_t)(dist_from_start / track_length * n_track_segments);
            if (segment >= n_track_segments)
                segment = n_track_segments - 1;
            //------------------------------------------------------------------
 
            for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
 
                // Electron gain in each GEM cascade (three draws always, as in the other methods)
                Double_t gain_gem1 = gRandom->Exp(15);
                Double_t gain_gem2 = gRandom->Exp(15);
                Double_t gain_gem3 = gRandom->Exp(15);
 
                if (gain_gem1 < 1.0)
                    gain_gem1 = 1.0;
                if (gain_gem2 < 1.0)
                    gain_gem2 = 1.0;
                if (gain_gem3 < 1.0)
                    gain_gem3 = 1.0;
 
                int total_gain = 0;
 
                if (zdist < level1) {
                    total_gain = 1.0;
                }
                if (zdist >= level1 && zdist < level2) {
                    total_gain = gain_gem3;
                }
                if (zdist >= level2 && zdist < level3) {
                    total_gain = gain_gem3 * gain_gem2;
                }
                if (zdist >= level3) {
                    total_gain = gain_gem3 * gain_gem2 * gain_gem1;
                }
 
                seg_local[segment] += total_gain;
 
                // ANN mode: electrons are not projected on the readout, strips are taken from the model
                if (useANN)
                    continue;
 
                // Standard method ---------------------------------------------
                // [PARALLEL-STD] Only the electron positions (exact mode) or the tasks (fast mode) are produced
                // here; strips are assigned later in parallel (see "Standard method: parallel strip assignment")
                // x-shift of the electon depends on the electron drift direction
                const Double_t x_center =
                    (ElectronDriftDirection == ForwardZAxisEDrift) ? (xcoll + xmean) : (xcoll - xmean);
 
                if (kExactGRandomSequence) {
                    // the same gRandom calls in the same order as in AddRealPointFull_ANNStripGenerator_MLP
                    for (int igain = 0; igain < total_gain; ++igain) {
                        exact_x.push_back(gRandom->Gaus(x_center, sigma));
                        exact_y.push_back(gRandom->Gaus(ycoll, sigma));
                    }
                } else {
                    // positions are generated later in parallel, here - only the tasks
                    for (Long64_t done = 0; done < total_gain; done += kElectronsPerTask) {
                        aval_tasks.push_back({x_center, ycoll, sigma, n_aval_electrons + done,
                                              std::min<Long64_t>(kElectronsPerTask, total_gain - done)});
                    }
                    n_aval_electrons += total_gain;
                }
            }
        }
    }
 
    // [PARALLEL-STD] Standard method: parallel strip assignment ------------
    // Each task counts its electrons per strip (compact local counters, no per-electron storage).
    // The counters of the tasks are summed in the task order and added to the clusters with one AddStrip()
    // per strip (as in AddRealPointFull_ANNStripGenerator_MLP): the signals are integer numbers of electrons,
    // so the clusters are bitwise the same as with one AddStrip() per electron and do not depend on the threads.
    if (!useANN) {
 
        // Number of avalanche electrons of this point
        const size_t n_electrons = kExactGRandomSequence ? exact_x.size() : static_cast<size_t>(n_aval_electrons);
 
        if (n_electrons > 0) {
 
            const Bool_t run_parallel = static_cast<Long64_t>(n_electrons) >= kMinElectronsForParallel;
            const size_t chunk_size = static_cast<size_t>(kElectronsPerTask);
 
            // Layer invariant (computed once instead of for each electron)
            std::vector<char> layer_vertical(n_strip_layers);
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
                layer_vertical[ilayer] = TMath::Abs(StripLayers[ilayer].GetAngleDeg()) < 1.0;
 
            // Counters of one task: in each layer the strips [lo, lo + cnt.size()) (the range grows if needed)
            struct TaskCounts
            {
                std::vector<Int_t> lo;
                std::vector<std::vector<Int_t>> cnt;
            };
            const size_t n_tasks =
                kExactGRandomSequence ? (n_electrons + chunk_size - 1) / chunk_size : aval_tasks.size();
            std::vector<TaskCounts> task_counts(n_tasks);
 
            // One electron -> +1 on its strip in each layer (the same conditions as in the old code)
            auto count_electron = [&](TaskCounts& tc, Double_t x_readout, Double_t y_readout) {
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                    // NOTE: these BmnGemStripLayer methods are called from several threads,
                    //       they must not modify the layer
                    // Condition: end electron position must be inside the layer
                    if (!StripLayers[ilayer].IsPointInsideStripLayer(x_readout, y_readout))
                        continue;
 
                    Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_readout, y_readout);
                    // AZ-100124 - account for strip width
                    Double_t dstrip = strip_pos - (int)strip_pos;
                    if (layer_vertical[ilayer]) {
                        // Vertical strips
                        if (dstrip >= 0.4 && dstrip < 0.6)
                            continue;
                    } else {
                        // Rotated strips
                        if (dstrip >= 0.3 && dstrip < 0.7)
                            continue;
                    }
                    // AZ
 
                    const Int_t strip = (Int_t)strip_pos;
                    std::vector<Int_t>& cnt = tc.cnt[ilayer];
                    Int_t& lo = tc.lo[ilayer];
                    if (cnt.empty()) {
                        lo = strip;
                        cnt.assign(1, 0);
                    } else if (strip < lo) {
                        cnt.insert(cnt.begin(), lo - strip, 0);
                        lo = strip;
                    } else if (strip >= lo + static_cast<Int_t>(cnt.size())) {
                        cnt.resize(strip - lo + 1, 0);
                    }
                    ++cnt[strip - lo];
                }
            };
 
            auto init_counts = [&](TaskCounts& tc) {
                tc.lo.assign(n_strip_layers, 0);
                tc.cnt.assign(n_strip_layers, std::vector<Int_t>());
            };
 
            if (kExactGRandomSequence) {
                // (a) Exact mode: positions are already generated (gRandom, sequentially) -> only geometry, by chunks
                auto process_chunk = [&](size_t ichunk) {
                    TaskCounts& tc = task_counts[ichunk];
                    init_counts(tc);
                    const size_t begin = ichunk * chunk_size;
                    const size_t end = std::min(n_electrons, begin + chunk_size);
                    for (size_t ie = begin; ie < end; ++ie)
                        count_electron(tc, exact_x[ie], exact_y[ie]);
                };
 
                if (run_parallel)
                    parallel_for(n_tasks, process_chunk);
                else
                    for (size_t ichunk = 0; ichunk < n_tasks; ++ichunk)
                        process_chunk(ichunk);
 
                // Electron points on the readout (generation data)
                if (FillGenerationData) {
                    XElectronPos.insert(XElectronPos.end(), exact_x.begin(), exact_x.end());
                    YElectronPos.insert(YElectronPos.end(), exact_y.begin(), exact_y.end());
                    ElectronSignal.resize(ElectronSignal.size() + n_electrons, 1.0);
                }
            } else {
                // (a) Fast mode: generation + geometry in one pass (one thread launch, one pass over memory)
 
                // base seed from gRandom -> reproducible with the seed of gRandom, independent of the threads
                const UInt_t base_seed = static_cast<UInt_t>(gRandom->Integer(0xFFFFFFFFu));
 
                // Electron points on the readout (generation data): each task writes its own range
                const size_t first_electron = XElectronPos.size();   // 0 after ResetElectronPointsAndClusters()
                if (FillGenerationData) {
                    XElectronPos.resize(first_electron + n_electrons);
                    YElectronPos.resize(first_electron + n_electrons);
                    ElectronSignal.resize(first_electron + n_electrons, 1.0);
                }
 
                auto process_task = [&](size_t itask) {
                    const AvalancheTask& t = aval_tasks[itask];
                    TaskCounts& tc = task_counts[itask];
                    init_counts(tc);
                    std::seed_seq seq{base_seed, static_cast<UInt_t>(itask)};
                    std::mt19937_64 rng(seq);
                    std::normal_distribution<Double_t> gaus(0.0, 1.0);   // mean + sigma*N(0,1), as TRandom::Gaus
 
                    for (Long64_t i = 0; i < t.n; ++i) {
                        const Double_t x_readout = t.xcenter + t.sigma * gaus(rng);
                        const Double_t y_readout = t.ycenter + t.sigma * gaus(rng);
 
                        if (FillGenerationData) {
                            const size_t ie = first_electron + static_cast<size_t>(t.offset + i);
                            XElectronPos[ie] = x_readout;
                            YElectronPos[ie] = y_readout;
                        }
 
                        count_electron(tc, x_readout, y_readout);
                    }
                };
 
                if (run_parallel)
                    parallel_for(n_tasks, process_task);
                else
                    for (size_t itask = 0; itask < n_tasks; ++itask)
                        process_task(itask);
            }
 
            // (b) Filling the clusters: the counters of all tasks are summed (in the task order),
            //     one AddStrip() per strip; layers are independent -> in parallel
            auto fill_layer = [&](size_t ilayer) {
                StripCluster& cluster = cluster_layers[ilayer];
                const Int_t n_strips_in_layer = std::max(0, StripLayers[ilayer].GetNStrips());
                std::vector<Double_t> strip_acc(n_strips_in_layer, 0.0);
                Int_t acc_min = std::numeric_limits<Int_t>::max();
                Int_t acc_max = -1;
 
                for (size_t itask = 0; itask < n_tasks; ++itask) {
                    const Int_t lo = task_counts[itask].lo[ilayer];
                    const std::vector<Int_t>& cnt = task_counts[itask].cnt[ilayer];
                    for (size_t i = 0; i < cnt.size(); ++i) {
                        if (cnt[i] == 0)
                            continue;
                        const Int_t strip = lo + static_cast<Int_t>(i);
                        if (strip >= 0 && strip < n_strips_in_layer) {
                            strip_acc[strip] += cnt[i];
                            if (strip < acc_min)
                                acc_min = strip;
                            if (strip > acc_max)
                                acc_max = strip;
                        } else {
                            // out of the accumulator range: directly (the same sum as cnt[i] x AddStrip(strip, 1.0))
                            cluster.AddStrip(strip, static_cast<Double_t>(cnt[i]));
                        }
                    }
                }
 
                for (Int_t istrip = acc_min; istrip <= acc_max; ++istrip) {
                    if (strip_acc[istrip] > 0.0)
                        cluster.AddStrip(istrip, strip_acc[istrip]);
                }
            };
 
            if (run_parallel)
                parallel_for(static_cast<size_t>(n_strip_layers), fill_layer);
            else
                for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer)
                    fill_layer(ilayer);
        }
    }
    //----------------------------------------------------------------------
 
    // ANN generation (MLP, own inference) ---------------------------------
    if (useANN) {
 
        // Segments for the model: external or from the avalanche of this point
        const vector<Int_t>& seg_used = use_external_segs ? segs : seg_local;
 
        // Normalization (px,py,pz), where r_vec must be equal to 1
        Double_t px_norm = px / r_vec;
        Double_t py_norm = py / r_vec;
        Double_t pz_norm = pz / r_vec;
 
        // Inversion for z-reverse direction particles: the same track segment passed from (x_out, y_out)
        if (pz < 0.0) {
            px_norm = -px_norm;
            py_norm = -py_norm;
            pz_norm = -pz_norm;
        }
 
        Double_t drift_direction = 0;   // 0 is ForwardZAxisEDrift
        if (ElectronDriftDirection == BackwardZAxisEDrift)
            drift_direction = 1;
 
        // Magnetic field: given explicitly or from the current medium configuration
        const Double_t bfield_used = (bfield < 0.0) ? BmnGemStripMedium::GetInstance().BField : bfield;
 
        // Normalize seg. vector: [0, max] -> [0, 1] (max over the whole vector)
        Double_t max_segment_val = *(std::max_element(seg_used.begin(), seg_used.end()));
 
        // Total charge in segment vector
        Double_t total_charge = 0.0;
        for (size_t i = 0; i < seg_used.size(); ++i)
            total_charge += seg_used[i];
 
        // Window center point (x_c, y_c): computed above with GetANNWindowCenterPoint()
 
        // Common part of the model input: [px, py, pz, drift, bfield, s0..s{n_segs_features-1}]
        float common_input[n_common_inputs];
        common_input[0] = static_cast<float>(px_norm);
        common_input[1] = static_cast<float>(py_norm);
        common_input[2] = static_cast<float>(pz_norm);
        common_input[3] = static_cast<float>(drift_direction);
        common_input[4] = static_cast<float>(bfield_used);
        for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
            common_input[5 + iseg] =
                (max_segment_val > 0.0) ? static_cast<float>(seg_used[iseg] / max_segment_val) : 0.0f;
 
        if (total_charge > 0.0) {
            // Strip model: a window of n_ann_strips strips around the window center strip in each layer
            const Int_t center_strip_index = n_ann_strips / 2;   // center strip of the window
 
            // Layers with the MC point inside: their part of the input [strip_layer, strip_order, center_frac]
            std::vector<Int_t> ann_layers;          // layer index
            std::vector<Int_t> ann_center_strips;   // MC strip (window center)
            std::vector<Int_t> ann_n_strips;        // number of strips in the layer
            std::vector<float> ann_layer_inputs;    // [n_ann_layers][n_layer_inputs]
            ann_layers.reserve(n_strip_layers);
            ann_center_strips.reserve(n_strip_layers);
            ann_n_strips.reserve(n_strip_layers);
            ann_layer_inputs.reserve(n_strip_layers * n_layer_inputs);
 
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
                // The window is built around the center point (middle of the drift gap): it must be inside this layer
                if (!StripLayers[ilayer].IsPointInsideStripLayer(x_c, y_c))
                    continue;
 
                Int_t n_strips_in_layer = StripLayers[ilayer].GetNStrips();
                if (n_strips_in_layer <= 0)
                    continue;
 
                // Center point position in strips: integer part - window center strip, fractional part - center_frac
                Double_t center_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_c, y_c);
                Int_t center_strip = (Int_t)TMath::Floor(center_pos);
                Double_t center_frac = center_pos - center_strip;
 
                Double_t strip_layer = static_cast<Double_t>(StripLayers[ilayer].GetType());   // 0 - lower, 1 - upper
                Double_t strip_order =
                    static_cast<Double_t>(StripLayers[ilayer].GetStripNumberingOrder());   // 0 - LeftToRight
 
                ann_layers.push_back(ilayer);
                ann_center_strips.push_back(center_strip);
                ann_n_strips.push_back(n_strips_in_layer);
                ann_layer_inputs.push_back(static_cast<float>(strip_layer));
                ann_layer_inputs.push_back(static_cast<float>(strip_order));
                ann_layer_inputs.push_back(static_cast<float>(center_frac));
            }
 
            // Strip signals for all these layers in one call: n_ann_strips values in [0, 1] (max ~ 1) per layer
            Int_t n_ann_layers = ann_layers.size();
            std::vector<float> ann_strip_values(static_cast<size_t>(n_ann_layers) * n_ann_strips);
            if (n_ann_layers > 0
                && !mlp_model->PredictSplit(common_input, n_common_inputs, ann_layer_inputs.data(), n_ann_layers,
                                            ann_strip_values.data()))
            {
                LOG(error) << "BmnGemStripModule::" << __func__ << "(): MLP prediction failed";
                n_ann_layers = 0;
            }
 
            for (Int_t iann = 0; iann < n_ann_layers; ++iann) {
                const Int_t ilayer = ann_layers[iann];
                const Int_t center_strip = ann_center_strips[iann];
                const Int_t n_strips_in_layer = ann_n_strips[iann];
                const float* strip_values = ann_strip_values.data() + static_cast<size_t>(iann) * n_ann_strips;
 
                // Sum of the model output over all strips of the window
                Double_t sum_norm = 0.0;
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip)
                    sum_norm += strip_values[istrip];
 
                // Scale coeff for total signal distribution (separate for lower and upper layers)
                const Double_t ann_correction_lower = 0.8;   // Set to 1.0 to disable correction
                const Double_t ann_correction_upper = 0.6;   // Set to 1.0 to disable correction
                const Double_t ann_correction =
                    (StripLayers[ilayer].GetType() == LowerStripLayer) ? ann_correction_lower : ann_correction_upper;
 
                Double_t signal_scale = 0.0;
                if (sum_norm > 0.0)
                    signal_scale = (total_charge / sum_norm) * ann_correction;
 
                Double_t threshold = 0.0001 * total_charge;
 
                // Strips are written directly: window index -> strip number
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip) {
                    Double_t content = strip_values[istrip] * signal_scale;
                    if (content <= threshold)
                        continue;
 
                    Int_t strip_num = center_strip - center_strip_index + istrip;
                    if (strip_num < 0 || strip_num >= n_strips_in_layer)
                        continue;   // the window goes out of the layer
 
                    cluster_layers[ilayer].AddStrip(strip_num, content);
                }
            }
        }
    }
    //----------------------------------------------------------------------
 
    // Segmented charge distribution (generation data) ------------------------
    // (seg_used above is not used after this point)
    if (FillGenerationData) {
        if (use_external_segs)
            SegmentedChargeDistribution = std::move(segs);
        else
            SegmentedChargeDistribution = std::move(seg_local);
    }
    //----------------------------------------------------------------------
 
    // Calculate cluster parameters --------------------------------------------
    vector<Double_t> cluster_mean_position(n_strip_layers, 0.0);
    vector<Double_t> cluster_total_signal(n_strip_layers, 0.0);
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
        if (cluster_layers[ilayer].GetClusterSize() > 0) {
 
            for (int i = 0; i < cluster_layers[ilayer].GetClusterSize(); ++i) {
                cluster_mean_position[ilayer] += (cluster_layers[ilayer].Strips[i] + 0.5)
                                                 * cluster_layers[ilayer].Signals[i];   // as sum of all positions
                cluster_total_signal[ilayer] += cluster_layers[ilayer].Signals[i];
            }
            cluster_mean_position[ilayer] /= cluster_total_signal[ilayer];
        }
 
        cluster_layers[ilayer].MeanPosition = cluster_mean_position[ilayer];   // mean cluster position
        cluster_layers[ilayer].TotalSignal = cluster_total_signal[ilayer];     // total signal of the cluster
    }
 
    if (ElectronDriftDirection == ForwardZAxisEDrift) {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        }
    } else {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        cluster_layers[ilayer].PositionResidual =
            cluster_layers[ilayer].MeanPosition - cluster_layers[ilayer].OriginPosition;
        cluster_layers[ilayer].SetType(StripLayers[ilayer].GetType());
    }
    //----------------------------------------------------------------------
 
    // Added clusters (generation data) ----------------------------------------
    if (FillGenerationData) {
        if (AddedClusters.size() == 0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                AddedClusters.push_back(vector<StripCluster>());
            }
        }
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
        }
    }
    //------------------------------------------------------------------------------
 
    // Add the correct clusters to the strip layers -------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
        }
    }
    //----------------------------------------------------------------------
 
    // Fill strip matches ---------------------------------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                    refID);
        }
    }
    //----------------------------------------------------------------------
 
    RealPointsX.push_back(x);
    RealPointsY.push_back(y);
    RealPointsXout.push_back(x_out);
    RealPointsYout.push_back(y_out);
    RealPointsMC.push_back(refID);
    RealPointsXcenter.push_back(x_c);   // ANN window center point
    RealPointsYcenter.push_back(y_c);
 
    return true;
}
//------------------------------------------------------------------------------
 
// The same as AddRealPointFull_ANNStripGenerator_MLP (ANN part: identical), but the standard method does not simulate
// the avalanche electrons one by one. All electrons of one ionizing collision (sum of the gains of its primary
// electrons) start from the same point and are smeared with the same Gaussian (Lorentz shift xmean, sigma).
// The strip coordinate of an electron is a linear function of (x, y), so in each strip layer it is Gaussian too,
// and the number of electrons in each strip is drawn from the multinomial distribution (a chain of binomials)
// with the exact strip probabilities (Gaussian integrals over the strip without its dead band).
// Statistically equivalent to the per-electron simulation; the layers are drawn independently (the correlation
// "the same electron in both layers" is lost: its effect on the cluster position is ~ sigma / sqrt(N electrons)).
// Collisions near a layer border (or at negative strip coordinates) are simulated per electron, as before.
// Random numbers: collisions, gains and segments - from gRandom in the same order as in the other strip methods;
// strips - from a std::mt19937_64 seeded from gRandom after the collisions (reproducible with the seed of gRandom).
// The electron positions (XElectronPos, YElectronPos) are not simulated and not filled.
Bool_t BmnGemStripModule::AddRealPointFull_ANNStripGenerator_MLP_FastStand(Double_t x,
                                                                           Double_t y,
                                                                           Double_t z,
                                                                           Double_t px,
                                                                           Double_t py,
                                                                           Double_t pz,
                                                                           Double_t signal,
                                                                           Int_t refID,
                                                                           Int_t pdg_code,
                                                                           Double_t bfield,
                                                                           Bool_t onlyStandard,
                                                                           vector<Int_t> segs)
{
    const Double_t maxThetaAngle = kANNMaxThetaAngle;   // max |theta| (deg) of the track for the ANN
 
    // Model format: [px, py, pz, drift, bfield, s0..s63, strip_layer, strip_order, center_frac] -> window of N
    // strips (N odd, from the model) centered at the middle of the drift gap (GetANNWindowCenterPoint)
    const Int_t n_track_segments = 64;                   // track segmentation (segmented charge distribution)
    const Int_t n_segs_features = n_track_segments;      // segments passed to the model
    const Int_t n_common_inputs = 5 + n_segs_features;   // the same for all strip layers of the point
    const Int_t n_layer_inputs = 3;                      // [strip_layer, strip_order, center_frac]
    Int_t n_ann_strips = 0;   // model output: strip window size (odd), taken from the model below
 
    // Test outputs are cleared always (cheap if they are empty)
    ResetElectronPointsAndClusters();
 
    // MLP model: shared between modules, checked once for each model ---------
    std::shared_ptr<const BmnGemStripMLPModel> mlp_holder;
    const BmnGemStripMLPModel* mlp_model = nullptr;
    if (!onlyStandard) {
        mlp_holder = BmnGemStripMLPModel::GetShared();
        mlp_model = mlp_holder.get();
 
        static std::mutex checked_mutex;
        static std::map<const BmnGemStripMLPModel*, std::pair<std::shared_ptr<const BmnGemStripMLPModel>, std::string>>
            checked;   // model -> (holder, why it is not suitable: empty if suitable)
 
        std::lock_guard<std::mutex> lock(checked_mutex);
        if (!mlp_model) {
            LOG(warn) << "BmnGemStripModule::" << __func__
                      << "(): MLP model is not available, the standard method is used";
        } else {
            auto it = checked.find(mlp_model);
            if (it == checked.end()) {
                std::vector<std::string> expected = {"px", "py", "pz", "drift", "bfield"};
                for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
                    expected.push_back("s" + std::to_string(iseg));
                expected.insert(expected.end(), {"strip_layer", "strip_order", "center_frac"});
 
                std::string why;
                Bool_t ok = mlp_model->CheckFeatureNames(expected, &why);
                if (ok) {
                    why = CheckANNStripWindowSize(mlp_model->GetOutputDim());
                    ok = why.empty();
                }
                if (ok) {
                    why.clear();
                    LOG(info) << "BmnGemStripModule::" << __func__ << "(): MLP model " << mlp_model->GetFilePath()
                              << " is used (ANN for |theta| < " << maxThetaAngle << " deg, window "
                              << mlp_model->GetOutputDim() << " strips)";
                } else if (why.empty())
                    why = "unknown reason";
                it = checked.emplace(mlp_model, std::make_pair(mlp_holder, why)).first;
            }
            if (!it->second.second.empty()) {
                LOG(warn) << "BmnGemStripModule::" << __func__ << "(): the model " << mlp_model->GetFilePath()
                          << " is not suitable (" << it->second.second << "), the standard method is used";
                mlp_model = nullptr;
            }
        }
    }
    if (mlp_model)
        n_ann_strips = mlp_model->GetOutputDim();   // window size of the model (checked above: odd, >= 3)
    //--------------------------------------------------------------------------
 
    // Condition: a start point is inside the module and outside its dead zone (in strip layers)
    if (!IsPointInsideModule(x, y)) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
 
    // Particle: pdg_code = 0 -> from the MC track of the point refID (as in AddRealPointFull) --
    Int_t pdg = pdg_code;
    if (pdg == 0) {
        FairRunSim* run = FairRunSim::Instance();
        BmnGemStripDigitizer* digi = run ? (BmnGemStripDigitizer*)run->GetTask("GemDigitizer") : nullptr;
        TClonesArray* points = digi ? digi->GetMCPoints() : nullptr;
        TClonesArray* mctracks = digi ? digi->GetMCTracks() : nullptr;
 
        CbmMCTrack* mcTr = nullptr;
        if (points && mctracks && refID >= 0 && refID < points->GetEntriesFast()) {
            CbmStsPoint* point = (CbmStsPoint*)points->UncheckedAt(refID);
            if (point && point->GetTrackID() >= 0 && point->GetTrackID() < mctracks->GetEntriesFast())
                mcTr = (CbmMCTrack*)mctracks->UncheckedAt(point->GetTrackID());
        }
 
        if (mcTr) {
            pdg = mcTr->GetPdgCode();
        } else {
            static std::atomic<bool> no_mc_reported{false};
            if (!no_mc_reported.exchange(true))
                LOG(warn) << "BmnGemStripModule::" << __func__
                          << "(): pdg_code = 0 and the MC track is not available (no GemDigitizer or MC point), "
                             "the He4 mass is used (this warning is printed once)";
        }
    }
 
    TVector3 mom(px, py, pz);
    TParticlePDG* part = (pdg != 0) ? TDatabasePDG::Instance()->GetParticle(pdg) : nullptr;
 
    Double_t mnoc = 1.0, mass = 3.7283999;   // He4
    if (part)
        mass = part->Mass();
    if (mass > 0.0001) {
        TLorentzVector lorv;
        lorv.SetVectM(mom, mass);
        Double_t beta = lorv.Beta();
        Double_t gamma = lorv.Gamma();
        Double_t charge = (part) ? part->Charge() / 3 : 1;
        mnoc = GetNumberOfClusters(beta, gamma, charge, 1.787, 12.33);
        if (mnoc < 0.5)
            return kFALSE;
        if (MaxClustersPerCm > 0.0 && mnoc > MaxClustersPerCm)
            mnoc = MaxClustersPerCm;
    }
    //--------------------------------------------------------------------------
 
    // Distance from entry point (x,y,z) to exit point (along z-axis)
    Double_t z_distance_from_in_to_out;
 
    // Condition: track direction along z-axis
    if (pz > 0.0) {
        z_distance_from_in_to_out = ModuleThickness - (z - ZStartModulePosition);
    } else {
        z_distance_from_in_to_out = z - ZStartModulePosition;
    }
 
    // Condition: track distance along z-axis must be not zero
    if (z_distance_from_in_to_out <= 0.0) {
        LOG(warn) << "BmnGemStripModule::" << __func__ << "(): Point (" << x << " : " << y << " : " << z
                  << ") has a track with an incorrect length";
        return false;
    }
 
    // Scale coefficient of the track vector (the path from entry point to exit point)
    Double_t vector_coeff = fabs(z_distance_from_in_to_out / pz);
 
    // Components of the track vector
    Double_t x_vec_from_in_to_out = vector_coeff * px;
    Double_t y_vec_from_in_to_out = vector_coeff * py;
    Double_t z_vec_from_in_to_out = vector_coeff * pz;
 
    // Exit point coordinates (x_out, y_out, z_out)
    Double_t x_out = x + x_vec_from_in_to_out;
    Double_t y_out = y + y_vec_from_in_to_out;
    Double_t z_out = z + z_vec_from_in_to_out;
 
    // Check if the exit point is outside the module then calculate new x_out, y_out, z_out values
    // Condition: x-coordinate of the exit point is inside the module
    if (x_out < XMinModule || x_out > XMaxModule) {
        if (x_out < XMinModule)
            x_out = XMinModule;
        if (x_out > XMaxModule)
            x_out = XMaxModule;
 
        x_vec_from_in_to_out = x_out - x;
        vector_coeff = x_vec_from_in_to_out / px;
 
        y_vec_from_in_to_out = vector_coeff * py;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        y_out = y + y_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    // Condition: y-coordinate of the exit point is inside the module
    if (y_out < YMinModule || y_out > YMaxModule) {
        if (y_out < YMinModule)
            y_out = YMinModule;
        if (y_out > YMaxModule)
            y_out = YMaxModule;
 
        y_vec_from_in_to_out = y_out - y;
        vector_coeff = y_vec_from_in_to_out / py;
 
        x_vec_from_in_to_out = vector_coeff * px;
        z_vec_from_in_to_out = vector_coeff * pz;
 
        x_out = x + x_vec_from_in_to_out;
        z_out = z + z_vec_from_in_to_out;
    }
    //--------------------------------------------------------------------------
 
    // Center point of the ANN strip window: the track at the middle of the drift gap (computed always, also with
    // onlyStandard: it is stored with the point, see GetRealPointXcenter(), for the training samples)
    Double_t x_c = x, y_c = y;
    GetANNWindowCenterPoint(x, y, z, x_out, y_out, z_out, x_c, y_c);
 
    // Common track length
    Double_t track_length =
        std::sqrt((x_vec_from_in_to_out) * (x_vec_from_in_to_out) + (y_vec_from_in_to_out) * (y_vec_from_in_to_out)
                  + (z_vec_from_in_to_out) * (z_vec_from_in_to_out));
 
    // Which method is used for this point: standard (electrons -> strips) or ANN (model -> strips)
    // (without a suitable model or with onlyStandard the standard method is used)
    Double_t r_vec = TMath::Sqrt(px * px + py * py + pz * pz);
    Double_t theta_deg_abs = TMath::ACos(TMath::Abs(pz) / r_vec) * TMath::RadToDeg();
    const Bool_t useANN = (theta_deg_abs < maxThetaAngle) && mlp_model;
 
    // ANN with external segments: the avalanche is not needed
    const Bool_t use_external_segs = useANN && !segs.empty();
    if (use_external_segs && segs.size() < static_cast<size_t>(n_segs_features)) {
        LOG(error) << "BmnGemStripModule::" << __func__ << "(): segment vector size = " << segs.size()
                   << ", the model requires " << n_segs_features;
        return false;
    }
 
    // Each level - distance to the readout plane
    Double_t level1 = InductionGapThickness;
    Double_t level2 = InductionGapThickness + SecondTransferGapThickness;
    Double_t level3 = InductionGapThickness + SecondTransferGapThickness + FirstTransferGapThickness;
 
    // Mean electron shift along x-axis (under the influence of the Lorentz force)
    Double_t xmean;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3)
 
    // Sigma electron smearing
    Double_t sigma;   // the dependence fitted by polynomial: f(x) = (p0 + p1*x + p2*x^2 + p3*x^3+ p4*x^4 + p5*x^5)
 
    // Reading the medium configuration -------------------------------------
    Double_t p0_xmean = BmnGemStripMedium::GetInstance().p0_xmean;
    Double_t p1_xmean = BmnGemStripMedium::GetInstance().p1_xmean;
    Double_t p2_xmean = BmnGemStripMedium::GetInstance().p2_xmean;
    Double_t p3_xmean = BmnGemStripMedium::GetInstance().p3_xmean;
 
    Double_t p0_sigma = BmnGemStripMedium::GetInstance().p0_sigma;
    Double_t p1_sigma = BmnGemStripMedium::GetInstance().p1_sigma;
    Double_t p2_sigma = BmnGemStripMedium::GetInstance().p2_sigma;
    Double_t p3_sigma = BmnGemStripMedium::GetInstance().p3_sigma;
    Double_t p4_sigma = BmnGemStripMedium::GetInstance().p4_sigma;
    Double_t p5_sigma = BmnGemStripMedium::GetInstance().p5_sigma;
    //----------------------------------------------------------------------
 
    // Initialization of a strip cluster in each strip layer
    Int_t n_strip_layers = StripLayers.size();   // number of layers in the module
    vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
    // Segmented charge distribution of this point (always computed when the avalanche is simulated)
    vector<Int_t> seg_local(n_track_segments, 0);
 
    // [FAST-STAND] Standard method: one task per ionizing collision ----------------
    // All electrons of a collision start from the same point -> one Gaussian (mean x_center, y_center; sigma);
    // n - sum of the gains of the primary electrons of the collision (number of electrons on the readout)
    struct CollisionTask
    {
        Double_t xcenter;   // mean x on the readout (with the Lorentz shift)
        Double_t ycenter;   // mean y on the readout
        Double_t sigma;     // smearing
        Long64_t n;         // number of electrons
    };
    std::vector<CollisionTask> coll_tasks;
    //--------------------------------------------------------------------------
 
    if (!use_external_segs) {
 
        Double_t current_length = 0.0;         // traversed length (counter)
        Double_t current_length_ratio = 0.0;   // ratio of the traversed length to common track length (counter)
        Double_t current_step = 0.0;           // current distance between two collision points
 
        // Collection of collision points
        std::vector<CollPoint> collision_points;
 
        // Reading the MCD value (mean collision distance = mean free flight path) [cm]
        MCD = BmnGemStripMedium::GetInstance().MCD;
        MCD = 1 / mnoc;   // AZ-110124
 
        while (current_length < track_length) {
 
            current_step = gRandom->Exp(MCD);
            current_length += current_step;
 
            if (current_length > track_length)
                break;
 
            current_length_ratio = current_length / track_length;
 
            Double_t current_x = x + current_length_ratio * x_vec_from_in_to_out;
            Double_t current_y = y + current_length_ratio * y_vec_from_in_to_out;
            Double_t current_z = z + current_length_ratio * z_vec_from_in_to_out;
 
            collision_points.push_back(CollPoint(current_x, current_y, current_z));
        }
 
        // Electron avalanche producing for each collision point at the track
        // (in the ANN mode only the segmented charge distribution is filled here)
        for (size_t icoll = 0; icoll < collision_points.size(); ++icoll) {
 
            Double_t xcoll = collision_points[icoll].x;
            Double_t ycoll = collision_points[icoll].y;
            Double_t zcoll = collision_points[icoll].z;
 
            Double_t zdist;   // current z-distance to the readout
 
            // Find z-distance to the readout depending on the electron drift direction
            if (ElectronDriftDirection == ForwardZAxisEDrift) {
                zdist = (ModuleThickness + ZStartModulePosition) - zcoll;
            } else {
                zdist = zcoll - ZStartModulePosition;
            }
            //------------------------------------------------------------------
 
            // polynomials
            xmean = p0_xmean + p1_xmean * zdist + p2_xmean * zdist * zdist + p3_xmean * zdist * zdist * zdist;
            sigma = p0_sigma + p1_sigma * zdist + p2_sigma * zdist * zdist + p3_sigma * zdist * zdist * zdist
                    + p4_sigma * zdist * zdist * zdist * zdist + p5_sigma * zdist * zdist * zdist * zdist * zdist;
 
            // Number of electrons in the current collision
            Int_t n_electrons_cluster = gRandom->Landau(1.027, 0.11);
            if (n_electrons_cluster < 1)
                n_electrons_cluster = 1;   // min
            if (n_electrons_cluster > 6)
                n_electrons_cluster = 6;   // max
 
            // Filling charge distribution -------------------------------------
            Double_t dist_from_start = 0.0;
            if (pz > 0.0) {
                dist_from_start =
                    std::sqrt((xcoll - x) * (xcoll - x) + (ycoll - y) * (ycoll - y) + (zcoll - z) * (zcoll - z));
            } else {
                dist_from_start = std::sqrt((xcoll - x_out) * (xcoll - x_out) + (ycoll - y_out) * (ycoll - y_out)
                                            + (zcoll - z_out) * (zcoll - z_out));
            }
            Int_t segment = (Int_t)(dist_from_start / track_length * n_track_segments);
            if (segment >= n_track_segments)
                segment = n_track_segments - 1;
            //------------------------------------------------------------------
 
            Long64_t coll_electrons = 0;   // [FAST-STAND] electrons of this collision on the readout
 
            for (Int_t ielectron = 0; ielectron < n_electrons_cluster; ++ielectron) {
 
                // Electron gain in each GEM cascade (three draws always, as in the other methods)
                Double_t gain_gem1 = gRandom->Exp(15);
                Double_t gain_gem2 = gRandom->Exp(15);
                Double_t gain_gem3 = gRandom->Exp(15);
 
                if (gain_gem1 < 1.0)
                    gain_gem1 = 1.0;
                if (gain_gem2 < 1.0)
                    gain_gem2 = 1.0;
                if (gain_gem3 < 1.0)
                    gain_gem3 = 1.0;
 
                int total_gain = 0;
 
                if (zdist < level1) {
                    total_gain = 1.0;
                }
                if (zdist >= level1 && zdist < level2) {
                    total_gain = gain_gem3;
                }
                if (zdist >= level2 && zdist < level3) {
                    total_gain = gain_gem3 * gain_gem2;
                }
                if (zdist >= level3) {
                    total_gain = gain_gem3 * gain_gem2 * gain_gem1;
                }
 
                seg_local[segment] += total_gain;
 
                // ANN mode: electrons are not projected on the readout, strips are taken from the model
                if (useANN)
                    continue;
 
                // Standard method: [FAST-STAND] the electrons of the collision are only counted here
                coll_electrons += total_gain;
            }
 
            // [FAST-STAND] one task per collision (the strips are drawn after the collision loop)
            if (!useANN && coll_electrons > 0) {
                const Double_t x_center =
                    (ElectronDriftDirection == ForwardZAxisEDrift) ? (xcoll + xmean) : (xcoll - xmean);
                coll_tasks.push_back({x_center, ycoll, sigma, coll_electrons});
            }
        }
    }
 
    // [FAST-STAND] Standard method: electrons -> strips by the multinomial distribution ----------
    // For each collision and each layer:
    //   u(x, y) = ConvertPointToStripPosition(x, y) is linear: u = u0 + a*(x - x0) + b*(y - y0)
    //   (x, y) ~ N(x0, sigma) x N(y0, sigma)  ->  u ~ N(u0, sigma_u),  sigma_u = sigma * sqrt(a^2 + b^2)
    //   strip j (u >= 0) gets the electrons with u in [j, j + d1) or [j + d2, j + 1) (d1, d2 - dead band of the strip,
    //   as in the per-electron code: (0.4, 0.6) for the vertical strips, (0.3, 0.7) for the rotated ones)
    //   numbers of electrons per strip ~ Multinomial(n, p_j): a chain of binomials, the rest (dead bands, far tails)
    //   is lost as in the per-electron code.
    // Collisions with the +-kNSigma box crossing the layer border, or with negative strip coordinates, are simulated
    // per electron (the same conditions as in the per-electron code).
    if (!useANN && !coll_tasks.empty()) {
 
        constexpr Double_t kNSigma = 6.0;   // range of the Gaussian (lost probability ~ 2e-9)
 
        // Generator of the strips: seeded from gRandom (reproducible with the seed of gRandom)
        std::mt19937_64 rng(static_cast<UInt_t>(gRandom->Integer(0xFFFFFFFFu)));
        std::normal_distribution<Double_t> gaus01(0.0, 1.0);
        std::binomial_distribution<Long64_t> binom;
        using BinomParam = std::binomial_distribution<Long64_t>::param_type;
 
        // Standard normal CDF
        auto Phi = [](Double_t t) { return 0.5 * TMath::Erfc(-t / TMath::Sqrt2()); };
 
        // Layer invariants: dead band, linear coefficients of the strip coordinate, accumulators
        struct LayerSetup
        {
            Bool_t linear = kFALSE;        // u(x, y) is linear (checked), else - per electron always
            Double_t d1 = 0.4, d2 = 0.6;   // dead band [d1, d2) of a strip (fraction of the strip)
            Double_t a = 0.0, b = 0.0;     // du/dx, du/dy
            Double_t grad = 0.0;           // sqrt(a^2 + b^2)
            Int_t n_strips = 0;
            std::vector<Double_t> acc;   // electrons per strip [0, n_strips)
            Int_t acc_min = std::numeric_limits<Int_t>::max(), acc_max = -1;
        };
        std::vector<LayerSetup> ls(n_strip_layers);
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            LayerSetup& s = ls[ilayer];
            BmnGemStripLayer& layer = StripLayers[ilayer];
            if (TMath::Abs(layer.GetAngleDeg()) >= 1.0) {
                s.d1 = 0.3;   // rotated strips
                s.d2 = 0.7;
            }
            s.n_strips = std::max(0, layer.GetNStrips());
            s.acc.assign(s.n_strips, 0.0);
 
            // linear coefficients from the strip coordinate around the window center point, h = 1 mm
            const Double_t h = 0.1;
            const Double_t u00 = layer.ConvertPointToStripPosition(x_c, y_c);
            const Double_t u10 = layer.ConvertPointToStripPosition(x_c + h, y_c);
            const Double_t u01 = layer.ConvertPointToStripPosition(x_c, y_c + h);
            const Double_t u11 = layer.ConvertPointToStripPosition(x_c + h, y_c + h);
            s.a = (u10 - u00) / h;
            s.b = (u01 - u00) / h;
            s.grad = TMath::Sqrt(s.a * s.a + s.b * s.b);
            // linearity check (the cross term must be zero)
            s.linear = std::isfinite(s.grad) && TMath::Abs(u11 - u10 - u01 + u00) <= 1e-9 * (1.0 + TMath::Abs(u00));
        }
 
        // One electron -> its strip in the layer (the same conditions as in the per-electron code)
        auto add_electron = [&](Int_t ilayer, Double_t x_readout, Double_t y_readout) {
            BmnGemStripLayer& layer = StripLayers[ilayer];
            if (!layer.IsPointInsideStripLayer(x_readout, y_readout))
                return;
            const Double_t strip_pos = layer.ConvertPointToStripPosition(x_readout, y_readout);
            const Double_t dstrip = strip_pos - (int)strip_pos;
            if (dstrip >= ls[ilayer].d1 && dstrip < ls[ilayer].d2)
                return;
            const Int_t strip = (Int_t)strip_pos;
            LayerSetup& s = ls[ilayer];
            if (strip >= 0 && strip < s.n_strips) {
                s.acc[strip] += 1.0;
                s.acc_min = std::min(s.acc_min, strip);
                s.acc_max = std::max(s.acc_max, strip);
            } else {
                cluster_layers[ilayer].AddStrip(strip, 1.0);   // out of the accumulator range: directly
            }
        };
 
        // n electrons in the strip -> accumulator
        auto add_strip = [&](Int_t ilayer, Int_t strip, Long64_t n) {
            LayerSetup& s = ls[ilayer];
            if (strip >= 0 && strip < s.n_strips) {
                s.acc[strip] += static_cast<Double_t>(n);
                s.acc_min = std::min(s.acc_min, strip);
                s.acc_max = std::max(s.acc_max, strip);
            } else {
                cluster_layers[ilayer].AddStrip(strip, static_cast<Double_t>(n));
            }
        };
 
        for (const CollisionTask& t : coll_tasks) {
            const Double_t box = kNSigma * t.sigma;   // half-size of the +-kNSigma box in x and y
 
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                LayerSetup& s = ls[ilayer];
                BmnGemStripLayer& layer = StripLayers[ilayer];
 
                // Position of the box relative to the layer: 4 corners + center
                Int_t n_inside = 0;
                n_inside += layer.IsPointInsideStripLayer(t.xcenter, t.ycenter) ? 1 : 0;
                n_inside += layer.IsPointInsideStripLayer(t.xcenter - box, t.ycenter - box) ? 1 : 0;
                n_inside += layer.IsPointInsideStripLayer(t.xcenter + box, t.ycenter - box) ? 1 : 0;
                n_inside += layer.IsPointInsideStripLayer(t.xcenter - box, t.ycenter + box) ? 1 : 0;
                n_inside += layer.IsPointInsideStripLayer(t.xcenter + box, t.ycenter + box) ? 1 : 0;
                if (n_inside == 0)
                    continue;   // the whole Gaussian is out of the layer
 
                const Double_t u0 = layer.ConvertPointToStripPosition(t.xcenter, t.ycenter);
                const Double_t sigma_u = t.sigma * s.grad;
                const Double_t u_lo = u0 - kNSigma * sigma_u;
                const Double_t u_hi = u0 + kNSigma * sigma_u;
 
                const Bool_t analytic = s.linear && n_inside == 5 && u_lo >= 0.0 && sigma_u > 1e-12;
 
                if (!analytic) {
                    // per electron (border of the layer, negative strip coordinates, no smearing)
                    for (Long64_t ie = 0; ie < t.n; ++ie) {
                        const Double_t x_readout = t.xcenter + t.sigma * gaus01(rng);
                        const Double_t y_readout = t.ycenter + t.sigma * gaus01(rng);
                        add_electron(ilayer, x_readout, y_readout);
                    }
                    continue;
                }
 
                // Multinomial over the strips j = floor(u_lo) .. floor(u_hi)
                const Int_t j_lo = (Int_t)TMath::Floor(u_lo);
                const Int_t j_hi = (Int_t)TMath::Floor(u_hi);
                Long64_t n_left = t.n;
                Double_t mass_left = 1.0;   // probability of the categories not drawn yet (incl. lost electrons)
 
                for (Int_t j = j_lo; j <= j_hi && n_left > 0; ++j) {
                    const Double_t p = (Phi((j + s.d1 - u0) / sigma_u) - Phi((j - u0) / sigma_u))
                                       + (Phi((j + 1 - u0) / sigma_u) - Phi((j + s.d2 - u0) / sigma_u));
                    if (p <= 0.0)
                        continue;
                    if (mass_left <= 0.0)
                        break;
                    Double_t q = p / mass_left;
                    if (q > 1.0)
                        q = 1.0;
                    const Long64_t k = binom(rng, BinomParam(n_left, q));
                    if (k > 0)
                        add_strip(ilayer, j, k);
                    n_left -= k;
                    mass_left -= p;
                }
                // n_left electrons: in the dead bands or out of +-kNSigma - lost (as in the per-electron code)
            }
        }
 
        // Accumulators -> clusters (one AddStrip() per strip, ascending order)
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            const LayerSetup& s = ls[ilayer];
            for (Int_t istrip = s.acc_min; istrip <= s.acc_max; ++istrip) {
                if (s.acc[istrip] > 0.0)
                    cluster_layers[ilayer].AddStrip(istrip, s.acc[istrip]);
            }
        }
    }
    //----------------------------------------------------------------------
 
    // ANN generation (MLP, own inference) ---------------------------------
    if (useANN) {
 
        // Segments for the model: external or from the avalanche of this point
        const vector<Int_t>& seg_used = use_external_segs ? segs : seg_local;
 
        // Normalization (px,py,pz), where r_vec must be equal to 1
        Double_t px_norm = px / r_vec;
        Double_t py_norm = py / r_vec;
        Double_t pz_norm = pz / r_vec;
 
        // Inversion for z-reverse direction particles: the same track segment passed from (x_out, y_out)
        if (pz < 0.0) {
            px_norm = -px_norm;
            py_norm = -py_norm;
            pz_norm = -pz_norm;
        }
 
        Double_t drift_direction = 0;   // 0 is ForwardZAxisEDrift
        if (ElectronDriftDirection == BackwardZAxisEDrift)
            drift_direction = 1;
 
        // Magnetic field: given explicitly or from the current medium configuration
        const Double_t bfield_used = (bfield < 0.0) ? BmnGemStripMedium::GetInstance().BField : bfield;
 
        // Normalize seg. vector: [0, max] -> [0, 1] (max over the whole vector)
        Double_t max_segment_val = *(std::max_element(seg_used.begin(), seg_used.end()));
 
        // Total charge in segment vector
        Double_t total_charge = 0.0;
        for (size_t i = 0; i < seg_used.size(); ++i)
            total_charge += seg_used[i];
 
        // Window center point (x_c, y_c): computed above with GetANNWindowCenterPoint()
 
        // Common part of the model input: [px, py, pz, drift, bfield, s0..s{n_segs_features-1}]
        float common_input[n_common_inputs];
        common_input[0] = static_cast<float>(px_norm);
        common_input[1] = static_cast<float>(py_norm);
        common_input[2] = static_cast<float>(pz_norm);
        common_input[3] = static_cast<float>(drift_direction);
        common_input[4] = static_cast<float>(bfield_used);
        for (Int_t iseg = 0; iseg < n_segs_features; ++iseg)
            common_input[5 + iseg] =
                (max_segment_val > 0.0) ? static_cast<float>(seg_used[iseg] / max_segment_val) : 0.0f;
 
        if (total_charge > 0.0) {
            // Strip model: a window of n_ann_strips strips around the window center strip in each layer
            const Int_t center_strip_index = n_ann_strips / 2;   // center strip of the window
 
            // Layers with the MC point inside: their part of the input [strip_layer, strip_order, center_frac]
            std::vector<Int_t> ann_layers;          // layer index
            std::vector<Int_t> ann_center_strips;   // MC strip (window center)
            std::vector<Int_t> ann_n_strips;        // number of strips in the layer
            std::vector<float> ann_layer_inputs;    // [n_ann_layers][n_layer_inputs]
            ann_layers.reserve(n_strip_layers);
            ann_center_strips.reserve(n_strip_layers);
            ann_n_strips.reserve(n_strip_layers);
            ann_layer_inputs.reserve(n_strip_layers * n_layer_inputs);
 
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
                // The window is built around the center point (middle of the drift gap): it must be inside this layer
                if (!StripLayers[ilayer].IsPointInsideStripLayer(x_c, y_c))
                    continue;
 
                Int_t n_strips_in_layer = StripLayers[ilayer].GetNStrips();
                if (n_strips_in_layer <= 0)
                    continue;
 
                // Center point position in strips: integer part - window center strip, fractional part - center_frac
                Double_t center_pos = StripLayers[ilayer].ConvertPointToStripPosition(x_c, y_c);
                Int_t center_strip = (Int_t)TMath::Floor(center_pos);
                Double_t center_frac = center_pos - center_strip;
 
                Double_t strip_layer = static_cast<Double_t>(StripLayers[ilayer].GetType());   // 0 - lower, 1 - upper
                Double_t strip_order =
                    static_cast<Double_t>(StripLayers[ilayer].GetStripNumberingOrder());   // 0 - LeftToRight
 
                ann_layers.push_back(ilayer);
                ann_center_strips.push_back(center_strip);
                ann_n_strips.push_back(n_strips_in_layer);
                ann_layer_inputs.push_back(static_cast<float>(strip_layer));
                ann_layer_inputs.push_back(static_cast<float>(strip_order));
                ann_layer_inputs.push_back(static_cast<float>(center_frac));
            }
 
            // Strip signals for all these layers in one call: n_ann_strips values in [0, 1] (max ~ 1) per layer
            Int_t n_ann_layers = ann_layers.size();
            std::vector<float> ann_strip_values(static_cast<size_t>(n_ann_layers) * n_ann_strips);
            if (n_ann_layers > 0
                && !mlp_model->PredictSplit(common_input, n_common_inputs, ann_layer_inputs.data(), n_ann_layers,
                                            ann_strip_values.data()))
            {
                LOG(error) << "BmnGemStripModule::" << __func__ << "(): MLP prediction failed";
                n_ann_layers = 0;
            }
 
            for (Int_t iann = 0; iann < n_ann_layers; ++iann) {
                const Int_t ilayer = ann_layers[iann];
                const Int_t center_strip = ann_center_strips[iann];
                const Int_t n_strips_in_layer = ann_n_strips[iann];
                const float* strip_values = ann_strip_values.data() + static_cast<size_t>(iann) * n_ann_strips;
 
                // Sum of the model output over all strips of the window
                Double_t sum_norm = 0.0;
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip)
                    sum_norm += strip_values[istrip];
 
                // Scale coeff for total signal distribution (separate for lower and upper layers)
                const Double_t ann_correction_lower = 0.8;   // Set to 1.0 to disable correction
                const Double_t ann_correction_upper = 0.6;   // Set to 1.0 to disable correction
                const Double_t ann_correction =
                    (StripLayers[ilayer].GetType() == LowerStripLayer) ? ann_correction_lower : ann_correction_upper;
 
                Double_t signal_scale = 0.0;
                if (sum_norm > 0.0)
                    signal_scale = (total_charge / sum_norm) * ann_correction;
 
                Double_t threshold = 0.0001 * total_charge;
 
                // Strips are written directly: window index -> strip number
                for (Int_t istrip = 0; istrip < n_ann_strips; ++istrip) {
                    Double_t content = strip_values[istrip] * signal_scale;
                    if (content <= threshold)
                        continue;
 
                    Int_t strip_num = center_strip - center_strip_index + istrip;
                    if (strip_num < 0 || strip_num >= n_strips_in_layer)
                        continue;   // the window goes out of the layer
 
                    cluster_layers[ilayer].AddStrip(strip_num, content);
                }
            }
        }
    }
    //----------------------------------------------------------------------
 
    // Segmented charge distribution (generation data) ------------------------
    // (seg_used above is not used after this point)
    if (FillGenerationData) {
        if (use_external_segs)
            SegmentedChargeDistribution = std::move(segs);
        else
            SegmentedChargeDistribution = std::move(seg_local);
    }
    //----------------------------------------------------------------------
 
    // Calculate cluster parameters --------------------------------------------
    vector<Double_t> cluster_mean_position(n_strip_layers, 0.0);
    vector<Double_t> cluster_total_signal(n_strip_layers, 0.0);
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
        if (cluster_layers[ilayer].GetClusterSize() > 0) {
 
            for (int i = 0; i < cluster_layers[ilayer].GetClusterSize(); ++i) {
                cluster_mean_position[ilayer] += (cluster_layers[ilayer].Strips[i] + 0.5)
                                                 * cluster_layers[ilayer].Signals[i];   // as sum of all positions
                cluster_total_signal[ilayer] += cluster_layers[ilayer].Signals[i];
            }
            cluster_mean_position[ilayer] /= cluster_total_signal[ilayer];
        }
 
        cluster_layers[ilayer].MeanPosition = cluster_mean_position[ilayer];   // mean cluster position
        cluster_layers[ilayer].TotalSignal = cluster_total_signal[ilayer];     // total signal of the cluster
    }
 
    if (ElectronDriftDirection == ForwardZAxisEDrift) {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        }
    } else {
        if (pz > 0.0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x_out, y_out);
            }
        } else {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                cluster_layers[ilayer].OriginPosition = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            }
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        cluster_layers[ilayer].PositionResidual =
            cluster_layers[ilayer].MeanPosition - cluster_layers[ilayer].OriginPosition;
        cluster_layers[ilayer].SetType(StripLayers[ilayer].GetType());
    }
    //----------------------------------------------------------------------
 
    // Added clusters (generation data) ----------------------------------------
    if (FillGenerationData) {
        if (AddedClusters.size() == 0) {
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                AddedClusters.push_back(vector<StripCluster>());
            }
        }
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
        }
    }
    //------------------------------------------------------------------------------
 
    // Add the correct clusters to the strip layers -------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
        }
    }
    //----------------------------------------------------------------------
 
    // Fill strip matches ---------------------------------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
            Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
            Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
            StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                    refID);
        }
    }
    //----------------------------------------------------------------------
 
    RealPointsX.push_back(x);
    RealPointsY.push_back(y);
    RealPointsXout.push_back(x_out);
    RealPointsYout.push_back(y_out);
    RealPointsMC.push_back(refID);
    RealPointsXcenter.push_back(x_c);   // ANN window center point
    RealPointsYcenter.push_back(y_c);
 
    return true;
}
//------------------------------------------------------------------------------
 
// Add single point with Gaussian smearing
Bool_t BmnGemStripModule::AddRealPointFullOne(Double_t x,
                                              Double_t y,
                                              Double_t z,
                                              Double_t px,
                                              Double_t py,
                                              Double_t pz,
                                              Double_t signal,
                                              Int_t refID)
{
 
    if (IsPointInsideModule(x, y)) {
 
        if (signal <= 0.0)
            signal = 1e-16;
 
        Double_t radius = AvalancheRadius;
        if (radius <= 0.0)
            return false;
 
        Int_t cycle_cnt = 0;
        while (true) {
            radius = gRandom->Gaus(AvalancheRadius, 0.02);
            if (radius > AvalancheRadius / 2.0 && radius < AvalancheRadius * 2.0 && radius > 0.0)
                break;
            cycle_cnt++;
 
            if (cycle_cnt > 5) {
                radius = AvalancheRadius;
                break;
            }
        }
 
        // Make strip cluster in each strip layer -------------------------------
        Int_t n_strip_layers = StripLayers.size();
        vector<StripCluster> cluster_layers(n_strip_layers, StripCluster());
 
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            if (!StripLayers[ilayer].IsPointInsideStripLayer(x, y))
                continue;
            cluster_layers[ilayer] = MakeCluster(ilayer, x, y, signal, radius);
        }
        //----------------------------------------------------------------------
 
        // Testing ----------------------------------------------------------------------
        /*if(AddedClusters.size() == 0) {
            for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
                AddedClusters.push_back(vector<StripCluster>());
            }
        }
 
        for(Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            AddedClusters[ilayer].push_back(cluster_layers[ilayer]);
        }*/
        //------------------------------------------------------------------------------
 
        // Add the correct clusters to the strip layers -------------------------
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
                Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
                Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
                StripLayers[ilayer].AddStripSignal(strip_num, strip_signal);
            }
        }
        //----------------------------------------------------------------------
 
        // Fill strip matches ---------------------------------------------------
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
            for (size_t ielement = 0; ielement < cluster_layers[ilayer].Strips.size(); ++ielement) {
                Int_t strip_num = cluster_layers[ilayer].Strips.at(ielement);
                Double_t strip_signal = cluster_layers[ilayer].Signals.at(ielement);
                StripLayers[ilayer].AddLinkToStripMatch(strip_num, strip_signal / cluster_layers[ilayer].TotalSignal,
                                                        refID);
            }
        }
        //----------------------------------------------------------------------
 
        RealPointsX.push_back(x);
        RealPointsY.push_back(y);
        RealPointsXout.push_back(x);
        RealPointsYout.push_back(y);
        RealPointsMC.push_back(refID);
 
        return true;
    } else {
        LOG(warn) << "BmnGemStripModule::AddRealPointFullOne(): Point (" << x << " : " << y << " : " << z
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
}
 
// Add single point without smearing and avalanch effects
Bool_t BmnGemStripModule::AddRealPointSimple(Double_t x,
                                             Double_t y,
                                             Double_t z,
                                             Double_t px,
                                             Double_t py,
                                             Double_t pz,
                                             Double_t signal,
                                             Int_t refID)
{
 
    if (IsPointInsideModule(x, y)) {
 
        Int_t n_strip_layers = StripLayers.size();
 
        // Add a signal to the strips layers ------------------------------------
        for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
            if (!StripLayers[ilayer].IsPointInsideStripLayer(x, y))
                continue;
 
            Double_t strip_pos = StripLayers[ilayer].ConvertPointToStripPosition(x, y);
            StripLayers[ilayer].AddStripSignal((Int_t)strip_pos, signal);
 
            // strip match
            StripLayers[ilayer].AddLinkToStripMatch((Int_t)strip_pos, 1.0, refID);
        }
        //----------------------------------------------------------------------
 
        RealPointsX.push_back(x);
        RealPointsY.push_back(y);
        RealPointsXout.push_back(x);
        RealPointsYout.push_back(y);
        RealPointsMC.push_back(refID);
 
        return true;
    } else {
        LOG(warn) << "BmnGemStripModule::AddRealPointSimple(): Point (" << x << " : " << y
                  << ") is out of the readout plane or inside a dead zone";
        return false;
    }
}
 
StripCluster BmnGemStripModule::MakeCluster(Int_t layer_num,
                                            Double_t xcoord,
                                            Double_t ycoord,
                                            Double_t signal,
                                            Double_t radius)
{
 
    Double_t ClusterDistortion = 0.0;   // signal variation (0.1 is 10%)
 
    StripCluster cluster;
 
    if (radius <= 0.0)
        return cluster;
 
    Double_t Pitch = StripLayers[layer_num].GetPitch();
 
    Double_t RadiusInZones = radius / Pitch;
    Double_t Sigma = RadiusInZones / 3.33;
 
    TF1 gausF("gausF", "gaus", -4 * Sigma, 4 * Sigma);
    gausF.SetParameter(0, 1.0);     // constant (altitude)
    gausF.SetParameter(1, 0.0);     // mean (center position)
    gausF.SetParameter(2, Sigma);   // sigma
 
    Double_t SCluster = gausF.Integral(-4 * Sigma, 4 * Sigma);   // square of the one side distribution (more exactly)
 
    TRandom rand(0);
    Double_t var_level = ClusterDistortion;   // signal variation (0.1 is 10%)
 
    Int_t NStripsInLayer = StripLayers[layer_num].GetNStrips();
 
    if (NStripsInLayer == 0)
        return cluster;
 
    Double_t CenterZonePos = StripLayers[layer_num].ConvertPointToStripPosition(xcoord, ycoord);
 
    if (CenterZonePos < 0.0 || CenterZonePos >= NStripsInLayer)
        return cluster;
 
    gausF.SetParameter(1, CenterZonePos);
    Double_t total_signal = 0.0;
 
    Double_t LeftZonePos = CenterZonePos - RadiusInZones;
    Double_t RightZonePos = CenterZonePos + RadiusInZones;
    Double_t OutLeftBorder = 0.0;
    Double_t OutRightBorder = 0.0;
 
    if (LeftZonePos < 0.0) {
        OutLeftBorder = LeftZonePos;
        LeftZonePos = 0.0;
    }
    if (RightZonePos < 0.0) {
        OutRightBorder = RightZonePos;
        RightZonePos = 0.0;
    }
    if (LeftZonePos >= NStripsInLayer) {
        OutLeftBorder = LeftZonePos;
        LeftZonePos = NStripsInLayer - 0.001;
    }
    if (RightZonePos >= NStripsInLayer) {
        OutRightBorder = RightZonePos;
        RightZonePos = NStripsInLayer - 0.001;
    }
 
    Double_t h = 0.0;
    Double_t dist = 0.0;
 
    // avalanche is inside of the one strip
    if ((Int_t)LeftZonePos == (Int_t)RightZonePos) {
 
        Int_t NumCurrentZone = (Int_t)LeftZonePos;
 
        h = RightZonePos - LeftZonePos;
        if (h < 0.0)
            h = 0.0;
 
        Double_t xp = LeftZonePos + dist;
        Double_t S = gausF.Integral(xp, xp + h);
        Double_t Energy = (signal * Gain) * (S / SCluster);
        Energy += rand.Gaus(0.0, var_level * Energy);
        if (Energy < 0.0)
            Energy = 0.0;
 
        if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
            cluster.AddStrip(NumCurrentZone, Energy);
            total_signal += Energy;
        }
 
        dist += h;
 
    } else {
        // left border strip
        Int_t NumCurrentZone = (Int_t)LeftZonePos;
 
        h = ((Int_t)LeftZonePos + 1) - LeftZonePos;
        if (h < 0.0)
            h = 0.0;
 
        Double_t xp = LeftZonePos + dist;
        Double_t S = gausF.Integral(xp, xp + h);
        Double_t Energy = (signal * Gain) * (S / SCluster);
        Energy += rand.Gaus(0.0, var_level * Energy);
        if (Energy < 0.0)
            Energy = 0.0;
 
        if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
            cluster.AddStrip(NumCurrentZone, Energy);
            total_signal += Energy;
        }
 
        dist += h;
 
        // inner strips
        for (Int_t i = (Int_t)LeftZonePos + 1; i < (Int_t)RightZonePos; ++i) {
 
            NumCurrentZone = i;
 
            h = 1.0;
 
            xp = LeftZonePos + dist;
            S = gausF.Integral(xp, xp + h);
            Energy = (signal * Gain) * (S / SCluster);
            Energy += rand.Gaus(0.0, var_level * Energy);
            if (Energy < 0.0)
                Energy = 0.0;
 
            if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
                cluster.AddStrip(NumCurrentZone, Energy);
                total_signal += Energy;
            }
 
            dist += h;
        }
        // right border strip
        NumCurrentZone = (Int_t)RightZonePos;
 
        h = RightZonePos - (Int_t)RightZonePos;
        if (h < 0.0)
            h = 0.0;
 
        xp = LeftZonePos + dist;
        S = gausF.Integral(xp, xp + h);
        Energy = (signal * Gain) * (S / SCluster);
        Energy += rand.Gaus(0.0, var_level * Energy);
        if (Energy < 0.0)
            Energy = 0.0;
 
        if (NumCurrentZone >= 0 && NumCurrentZone < NStripsInLayer && Energy > 0.0) {
            cluster.AddStrip(NumCurrentZone, Energy);
            total_signal += Energy;
        }
 
        dist += h;
    }
 
    if (cluster.GetClusterSize() <= 0) {
        return cluster;
    }
 
    // find the mean value of the avalanche position (fitting by gaus function)
    Double_t mean_fit_pos = 0.0;
 
    Int_t NOutLeftBins = 0;
    Int_t NOutRightBins = 0;
    if (OutLeftBorder != 0.0) {
        NOutLeftBins = (Int_t)(fabs(OutLeftBorder)) + 1;
    }
    if (OutRightBorder != 0.0) {
        NOutRightBins = (Int_t)(OutRightBorder - RightZonePos) + 1;
    }
 
    Int_t begin_strip_num = cluster.Strips.at(0);
    Int_t last_strip_num = cluster.Strips.at(cluster.GetClusterSize() - 1);
    Int_t nstrips = last_strip_num - begin_strip_num + 1;
 
    TH1F hist("hist_for_fit", "hist_for_fit", nstrips + NOutLeftBins + NOutRightBins, begin_strip_num - NOutLeftBins,
              last_strip_num + 1 + NOutRightBins);
    Int_t hist_index = 0;
 
    for (Int_t i = 0; i < cluster.GetClusterSize(); ++i) {
        Double_t value = cluster.Signals.at(i);
        hist.SetBinContent(hist_index + 1 + NOutLeftBins, value);
        hist_index++;
    }
 
    // on the left border
    if (NOutLeftBins > 0.0) {
        Double_t first = OutLeftBorder;
        Double_t last = (Int_t)OutLeftBorder;
        Double_t S = gausF.Integral(first, last);
        Double_t Energy = (signal * Gain) * (S / SCluster);
        Energy += rand.Gaus(0.0, var_level * Energy);
        if (Energy < 0.0)
            Energy = 0.0;
        Double_t value = Energy;
        hist.SetBinContent(1, value);
        dist = 0.0;
 
        for (Int_t i = 1; i < NOutLeftBins; ++i) {
            h = 1.0;
            first = (Int_t)OutLeftBorder + dist;
            last = first + h;
            S = gausF.Integral(first, last);
            Energy = (signal * Gain) * (S / SCluster);
            Energy += rand.Gaus(0.0, var_level * Energy);
            if (Energy < 0.0)
                Energy = 0.0;
            value = Energy;
            hist.SetBinContent(1 + i, value);
            dist += h;
        }
    }
 
    // on the right border
    if (NOutRightBins > 0.0) {
        dist = 0.0;
        for (Int_t i = hist_index; i < hist_index + NOutRightBins - 1; ++i) {
            h = 1.0;
            Double_t first = NStripsInLayer + dist;
            Double_t last = first + h;
            Double_t S = gausF.Integral(first, last);
            Double_t Energy = (signal * Gain) * (S / SCluster);
            Energy += rand.Gaus(0.0, var_level * Energy);
            if (Energy < 0.0)
                Energy = 0.0;
            Double_t value = Energy;
            hist.SetBinContent(1 + i, value);
            dist += h;
        }
 
        Double_t first = (Int_t)OutRightBorder;
        Double_t last = first + (OutRightBorder - (Int_t)OutRightBorder);
        Double_t S = gausF.Integral(first, last);
        Double_t Energy = (signal * Gain) * (S / SCluster);
        Energy += rand.Gaus(0.0, var_level * Energy);
        if (Energy < 0.0)
            Energy = 0.0;
        Double_t value = Energy;
        hist.SetBinContent(hist_index + NOutRightBins, value);
    }
 
    TF1* gausFitFunction = 0;
    TString fit_params = "WQ0";
 
#ifdef DRAW_REAL_CLUSTER_HISTOGRAMS
    fit_params = "WQ";
#endif
 
    if (nstrips > 1) {
        hist.Fit("gaus", fit_params);   // Q - quit mode (without information on the screen); 0 - not draw
        gausFitFunction = hist.GetFunction("gaus");
        if (gausFitFunction) {
            mean_fit_pos = gausFitFunction->GetParameter(1);
        } else {
            mean_fit_pos = hist.GetMean();
        }
    } else {
        mean_fit_pos = hist.GetMean();
    }
 
    cluster.MeanPosition = mean_fit_pos;
    cluster.TotalSignal = total_signal;
    cluster.PositionResidual = mean_fit_pos - CenterZonePos;   // residual between mean and origin positions (lower
                                                               // cluster): residual = finded(current) - orig(average)
 
    cluster.IsCorrect = true;   // set correct status of the cluster
    return cluster;
}
 
void BmnGemStripModule::SetAvalancheRadius(Double_t aval_radius)
{
    if (aval_radius >= 0.0)
        AvalancheRadius = aval_radius;
    else {
        LOG(warn) << "BmnGemStripModule::SetAvalancheRadius(): incorrect value for AvalancheRadius! It'll be "
                     "set by default!";
        AvalancheRadius = 0.1;   // default
    }
}
 
Double_t BmnGemStripModule::GetAvalancheRadius()
{
    return AvalancheRadius;
}
 
void BmnGemStripModule::CalculateStripHitIntersectionPoints()
{
 
    ResetIntersectionPoints();
    ResetPseudoIntersections();
 
    Int_t n_strip_layers = StripLayers.size();
 
    // Find strip clusters and hits in the strip layers -------------------------
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        StripLayers[ilayer].FindClustersAndStripHits();
    }
    //--------------------------------------------------------------------------
 
    map<Int_t, Bool_t> UsageStatus_LowerClusters;
    map<Int_t, Bool_t> UsageStatus_UpperClusters;
 
    // Saving all unique id marks of clusters in a module
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (Int_t i_strip_hit = 0; i_strip_hit < StripLayers[ilayer].GetNStripHits(); ++i_strip_hit) {
            StripCluster cluster = StripLayers[ilayer].GetStripClusters()[i_strip_hit];
            if (cluster.GetType() == LowerStripLayer) {
                UsageStatus_LowerClusters[cluster.GetUniqueID()] = false;
            }
            if (cluster.GetType() == UpperStripLayer) {
                UsageStatus_UpperClusters[cluster.GetUniqueID()] = false;
            }
        }
    }
 
    for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
        for (Int_t jlayer = ilayer + 1; jlayer < n_strip_layers; ++jlayer) {
            // cout << "(i-layer : j-layer) = ( " << ilayer << " : " << jlayer << " )\n";
 
            if (StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType())
                continue;
 
            Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
            Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
 
            if (Abs(iAngleRad - jAngleRad) < 0.01)
                continue;
 
            Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
            Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
 
            for (Int_t i_strip_hit = 0; i_strip_hit < StripLayers[ilayer].GetNStripHits(); ++i_strip_hit) {
                for (Int_t j_strip_hit = 0; j_strip_hit < StripLayers[jlayer].GetNStripHits(); ++j_strip_hit) {
 
                    Double_t iStripHitPos = StripLayers[ilayer].GetStripHitPos(i_strip_hit);
                    Double_t jStripHitPos = StripLayers[jlayer].GetStripHitPos(j_strip_hit);
 
                    Double_t iStripHitError = StripLayers[ilayer].GetStripHitError(i_strip_hit);
                    Double_t jStripHitError = StripLayers[jlayer].GetStripHitError(j_strip_hit);
 
                    Double_t xcoord;
                    Double_t ycoord;
 
                    if (Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01) {
 
                        Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
                        Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
 
                        xcoord = (jB - iB) / (Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
                        ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
                    } else {
 
                        if (Abs(iAngleRad) < 0.01) {
 
                            if (StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight) {
                                xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint()
                                         + iStripHitPos * StripLayers[ilayer].GetPitch();
                            } else {
                                xcoord = StripLayers[ilayer].GetXRightStripBorderPoint()
                                         - iStripHitPos * StripLayers[ilayer].GetPitch();
                            }
 
                            Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
                            ycoord = Tan(j_PiOver2MinusAngleRad) * xcoord + jB;
                        }
 
                        if (Abs(jAngleRad) < 0.01) {
 
                            if (StripLayers[jlayer].GetStripNumberingOrder() == LeftToRight) {
                                xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint()
                                         + jStripHitPos * StripLayers[jlayer].GetPitch();
                            } else {
                                xcoord = StripLayers[jlayer].GetXRightStripBorderPoint()
                                         - jStripHitPos * StripLayers[jlayer].GetPitch();
                            }
 
                            Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
                            ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
                        }
                    }
 
                    // check: a point belongs to both strip layers together
                    if (StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord)
                        && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord))
                    {
 
                        IntersectionPointsX.push_back(xcoord);
                        IntersectionPointsY.push_back(ycoord);
 
                        // Intersection (x,y)-errors from the strip-errors ------
                        Double_t AngleIntersecRad = Abs(iAngleRad - jAngleRad);
 
                        // a rhomb zone is the intersection of the strip errors
                        Double_t RhombSide1 = (iStripHitError * StripLayers[ilayer].GetPitch()) / Sin(AngleIntersecRad);
                        Double_t RhombSide2 = (jStripHitError * StripLayers[jlayer].GetPitch()) / Sin(AngleIntersecRad);
 
                        // x-strip error
                        Double_t xerr1 = Sin(Abs(jAngleRad)) * RhombSide1;
                        Double_t xerr2 = Sin(Abs(iAngleRad)) * RhombSide2;
                        Double_t xcoord_err = xerr1 + xerr2;
 
                        // y-strip error
                        Double_t yerr1 = Cos(Abs(jAngleRad)) * RhombSide1;
                        Double_t yerr2 = Cos(Abs(iAngleRad)) * RhombSide2;
                        Double_t ycoord_err = yerr1 + yerr2;
 
                        IntersectionPointsXErrors.push_back(xcoord_err);
                        IntersectionPointsYErrors.push_back(ycoord_err);
                        //------------------------------------------------------
 
                        // intersection MC-matching -----------------------------
                        BmnMatch iStripMatch = StripLayers[ilayer].GetStripMatch((Int_t)iStripHitPos);
                        BmnMatch jStripMatch = StripLayers[jlayer].GetStripMatch((Int_t)jStripHitPos);
 
                        BmnMatch intersection_match;
                        intersection_match.AddLink(iStripMatch);
                        intersection_match.AddLink(jStripMatch);
 
                        IntersectionPointMatches.push_back(intersection_match);
                        //------------------------------------------------------
 
                        // intersection digit number matching -------------------
                        BmnMatch iStripDigitNumberMatch =
                            StripLayers[ilayer].GetStripDigitNumberMatch((Int_t)iStripHitPos);
                        BmnMatch jStripDigitNumberMatch =
                            StripLayers[jlayer].GetStripDigitNumberMatch((Int_t)jStripHitPos);
 
                        BmnMatch intersection_digit_num_match;
                        intersection_digit_num_match.AddLink(iStripDigitNumberMatch);
                        intersection_digit_num_match.AddLink(jStripDigitNumberMatch);
 
                        IntersectionPointDigitNumberMatches.push_back(intersection_digit_num_match);
                        //------------------------------------------------------
 
                        // Additional information about the intersection ------------
                        // cluster size (number strips) in both strip layers for each intersection
                        Int_t iLayerClusterSize = StripLayers[ilayer].GetStripHitClusterSize(i_strip_hit);
                        Int_t jLayerClusterSize = StripLayers[jlayer].GetStripHitClusterSize(j_strip_hit);
 
                        // strip position in both strip layers for each intersection
                        Double_t iLayerStripPosition = StripLayers[ilayer].ConvertPointToStripPosition(xcoord, ycoord);
                        Double_t jLayerStripPosition = StripLayers[jlayer].ConvertPointToStripPosition(xcoord, ycoord);
 
                        if (StripLayers[ilayer].GetType() == LowerStripLayer) {
                            IntersectionPoints_LowerLayerClusterSize.push_back(iLayerClusterSize);
                            IntersectionPoints_UpperLayerClusterSize.push_back(jLayerClusterSize);
                            IntersectionPoints_LowerLayerStripPosition.push_back(iLayerStripPosition);
                            IntersectionPoints_UpperLayerStripPosition.push_back(jLayerStripPosition);
                            IntersectionPoints_LowerLayerStripTotalSignal.push_back(
                                StripLayers[ilayer].GetStripHitTotalSignal(i_strip_hit));
                            IntersectionPoints_UpperLayerStripTotalSignal.push_back(
                                StripLayers[jlayer].GetStripHitTotalSignal(j_strip_hit));
                            LowerClusters.push_back(StripLayers[ilayer].GetStripClusters().at(i_strip_hit));
                            UpperClusters.push_back(StripLayers[jlayer].GetStripClusters().at(j_strip_hit));
                            UsageStatus_LowerClusters
                                [StripLayers[ilayer].GetStripClusters()[i_strip_hit].GetUniqueID()] = true;
                            UsageStatus_UpperClusters
                                [StripLayers[jlayer].GetStripClusters()[j_strip_hit].GetUniqueID()] = true;
                        } else {
                            IntersectionPoints_LowerLayerClusterSize.push_back(jLayerClusterSize);
                            IntersectionPoints_UpperLayerClusterSize.push_back(iLayerClusterSize);
                            IntersectionPoints_LowerLayerStripPosition.push_back(jLayerStripPosition);
                            IntersectionPoints_UpperLayerStripPosition.push_back(iLayerStripPosition);
                            IntersectionPoints_LowerLayerStripTotalSignal.push_back(
                                StripLayers[jlayer].GetStripHitTotalSignal(j_strip_hit));
                            IntersectionPoints_UpperLayerStripTotalSignal.push_back(
                                StripLayers[ilayer].GetStripHitTotalSignal(i_strip_hit));
                            LowerClusters.push_back(StripLayers[jlayer].GetStripClusters().at(j_strip_hit));
                            UpperClusters.push_back(StripLayers[ilayer].GetStripClusters().at(i_strip_hit));
                            UsageStatus_LowerClusters
                                [StripLayers[jlayer].GetStripClusters()[j_strip_hit].GetUniqueID()] = true;
                            UsageStatus_UpperClusters
                                [StripLayers[ilayer].GetStripClusters()[i_strip_hit].GetUniqueID()] = true;
                        }
                        //------------------------------------------------------
                    }
                }
            }
        }
    }
    /*
        for(auto it : UsageStatus_UpperClusters) {
            if(it.second == false) cout << " upper_cluster[" << it.first << "] = " << it.second << "\n";
        }
    */
    for (auto it : UsageStatus_LowerClusters) {
 
        // if(it.second == false) cout << " lower_cluster[" << it.first << "] = " << it.second << "\n";
 
        if (it.second == false) {
 
            for (Int_t ilayer = 0; ilayer < n_strip_layers; ++ilayer) {
 
                if (StripLayers[ilayer].GetType() != LowerStripLayer)
                    continue;
 
                for (size_t icluster = 0; icluster < StripLayers[ilayer].GetStripClusters().size(); ++icluster) {
                    StripCluster cluster = StripLayers[ilayer].GetStripClusters()[icluster];
 
                    if ((int)cluster.GetUniqueID() != it.first)
                        continue;
 
                    Double_t xcoord;
                    Double_t ycoord;
 
                    if (StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight) {
                        xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint()
                                 + cluster.MeanPosition * StripLayers[ilayer].GetPitch();
                    } else {
                        xcoord = StripLayers[ilayer].GetXRightStripBorderPoint()
                                 - cluster.MeanPosition * StripLayers[ilayer].GetPitch();
                    }
 
                    ycoord = (StripLayers[ilayer].GetYMinLayer() + StripLayers[ilayer].GetYMaxLayer()) * 0.5;
 
                    Double_t xerr = cluster.Error * StripLayers[ilayer].GetPitch();
                    Double_t yerr = StripLayers[ilayer].GetYSize() * 0.5;
 
                    // cout << " low_cluster[" << icluster << "]:\n";
                    // cout << "  ID = " << cluster.GetUniqueID() << "\n";
                    // cout << "  xcoord = " << xcoord << "\n";
                    // cout << "   xmin = " << StripLayers[ilayer].GetXMinLayer() << ", xmax = " <<
                    // StripLayers[ilayer].GetXMaxLayer() << "\n"; cout << "  ycoord = " << ycoord << "\n"; cout << "
                    // ymin = " << StripLayers[ilayer].GetYMinLayer() << ", ymax = " <<
                    // StripLayers[ilayer].GetYMaxLayer() << "\n"; cout << "  xerr = " << xerr << "\n"; cout << "  yerr
                    // = " << yerr << "\n"; cout << "\n";
 
                    PseudoIntersectionsX.push_back(xcoord);
                    PseudoIntersectionsY.push_back(ycoord);
 
                    PseudoIntersectionsXErrors.push_back(xerr);
                    PseudoIntersectionsYErrors.push_back(yerr);
 
                    PseudoIntersections_LowerLayerClusterSize.push_back(cluster.Width);
                    PseudoIntersections_UpperLayerClusterSize.push_back(0);
 
                    PseudoIntersections_LowerLayerStripPosition.push_back(
                        StripLayers[ilayer].ConvertPointToStripPosition(xcoord, ycoord));
                    PseudoIntersections_UpperLayerStripPosition.push_back(0);
 
                    PseudoIntersections_LowerLayerStripTotalSignal.push_back(cluster.TotalSignal);
                    PseudoIntersections_UpperLayerStripTotalSignal.push_back(0);
 
                    PseudoIntersectionMatches.push_back(
                        StripLayers[ilayer].GetStripMatch((Int_t)cluster.MeanPosition));   // mc-match
                    PseudoIntersectionDigitNumberMatches.push_back(StripLayers[ilayer].GetStripDigitNumberMatch(
                        (Int_t)cluster.MeanPosition));   // digit number matches
 
                    LowerClusters_PseudoIntersections.push_back(cluster);
                    UpperClusters_PseudoIntersections.push_back(
                        StripCluster());   // empty cluster with default parameters
                }
            }
        }
    }
}
 
// need for a separated test (find x,y intersection coords from strip positions)
Bool_t BmnGemStripModule::SearchIntersectionPoint(Double_t& x,
                                                  Double_t& y,
                                                  Double_t strip_pos_layerA,
                                                  Double_t strip_pos_layerB,
                                                  Int_t layerA_index,
                                                  Int_t layerB_index)
{
 
    Int_t ilayer = layerA_index;
    Int_t jlayer = layerB_index;
 
    if (StripLayers[ilayer].GetType() == StripLayers[jlayer].GetType())
        return false;
 
    Double_t iAngleRad = StripLayers[ilayer].GetAngleRad();
    Double_t jAngleRad = StripLayers[jlayer].GetAngleRad();
 
    if (Abs(iAngleRad - jAngleRad) < 0.01)
        return false;
 
    Double_t i_PiOver2MinusAngleRad = PiOver2() - iAngleRad;
    Double_t j_PiOver2MinusAngleRad = PiOver2() - jAngleRad;
 
    Double_t iStripHitPos = strip_pos_layerA;
    Double_t jStripHitPos = strip_pos_layerB;
 
    Double_t xcoord;
    Double_t ycoord;
 
    if (Abs(iAngleRad) >= 0.01 && Abs(jAngleRad) >= 0.01) {
 
        Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
        Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
 
        xcoord = (jB - iB) / (Tan(i_PiOver2MinusAngleRad) - Tan(j_PiOver2MinusAngleRad));
        ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
    } else {
 
        if (Abs(iAngleRad) < 0.01) {
 
            if (StripLayers[ilayer].GetStripNumberingOrder() == LeftToRight) {
                xcoord = StripLayers[ilayer].GetXLeftStripBorderPoint() + iStripHitPos * StripLayers[ilayer].GetPitch();
            } else {
                xcoord =
                    StripLayers[ilayer].GetXRightStripBorderPoint() - iStripHitPos * StripLayers[ilayer].GetPitch();
            }
 
            Double_t jB = StripLayers[jlayer].CalculateStripEquationB(jStripHitPos);
            ycoord = Tan(j_PiOver2MinusAngleRad) * xcoord + jB;
        }
 
        if (Abs(jAngleRad) < 0.01) {
 
            if (StripLayers[jlayer].GetStripNumberingOrder() == LeftToRight) {
                xcoord = StripLayers[jlayer].GetXLeftStripBorderPoint() + jStripHitPos * StripLayers[jlayer].GetPitch();
            } else {
                xcoord =
                    StripLayers[jlayer].GetXRightStripBorderPoint() - jStripHitPos * StripLayers[jlayer].GetPitch();
            }
 
            Double_t iB = StripLayers[ilayer].CalculateStripEquationB(iStripHitPos);
            ycoord = Tan(i_PiOver2MinusAngleRad) * xcoord + iB;
        }
    }
 
    x = xcoord;
    y = ycoord;
 
    // check: a point belongs to both strip layers together
    if (StripLayers[ilayer].IsPointInsideStripLayer(xcoord, ycoord)
        && StripLayers[jlayer].IsPointInsideStripLayer(xcoord, ycoord))
    {
        return true;
    }
 
    return false;
}
 
void BmnGemStripModule::ResetIntersectionPoints()
{
    IntersectionPointsX.clear();
    IntersectionPointsY.clear();
    IntersectionPointsXErrors.clear();
    IntersectionPointsYErrors.clear();
    IntersectionPoints_LowerLayerClusterSize.clear();
    IntersectionPoints_UpperLayerClusterSize.clear();
    IntersectionPoints_LowerLayerStripPosition.clear();
    IntersectionPoints_UpperLayerStripPosition.clear();
    IntersectionPoints_LowerLayerStripTotalSignal.clear();
    IntersectionPoints_UpperLayerStripTotalSignal.clear();
    IntersectionPointMatches.clear();
    IntersectionPointDigitNumberMatches.clear();
    UpperClusters.clear();
    LowerClusters.clear();
}
 
void BmnGemStripModule::ResetPseudoIntersections()
{
    PseudoIntersectionsX.clear();
    PseudoIntersectionsY.clear();
    PseudoIntersectionsXErrors.clear();
    PseudoIntersectionsYErrors.clear();
    PseudoIntersections_LowerLayerClusterSize.clear();
    PseudoIntersections_UpperLayerClusterSize.clear();
    PseudoIntersections_LowerLayerStripPosition.clear();
    PseudoIntersections_UpperLayerStripPosition.clear();
    PseudoIntersections_LowerLayerStripTotalSignal.clear();
    PseudoIntersections_UpperLayerStripTotalSignal.clear();
    PseudoIntersectionMatches.clear();
    PseudoIntersectionDigitNumberMatches.clear();
    UpperClusters_PseudoIntersections.clear();
    LowerClusters_PseudoIntersections.clear();
}
 
void BmnGemStripModule::DefineModuleBorders()
{
 
    if (StripLayers.size() == 0)
        return;
 
    XMinModule = 1.0E+10;
    XMaxModule = -1.0E+10;
    YMinModule = 1.0E+10;
    YMaxModule = -1.0E+10;
 
    for (size_t ilayer = 0; ilayer < StripLayers.size(); ++ilayer) {
        if (StripLayers[ilayer].GetXMinLayer() < XMinModule)
            XMinModule = StripLayers[ilayer].GetXMinLayer();
        if (StripLayers[ilayer].GetXMaxLayer() > XMaxModule)
            XMaxModule = StripLayers[ilayer].GetXMaxLayer();
        if (StripLayers[ilayer].GetYMinLayer() < YMinModule)
            YMinModule = StripLayers[ilayer].GetYMinLayer();
        if (StripLayers[ilayer].GetYMaxLayer() > YMaxModule)
            YMaxModule = StripLayers[ilayer].GetYMaxLayer();
    }
    return;
}
 
// -------------------------------------------------------------------------
 
Double_t BmnGemStripModule::GetNumberOfClusters(Double_t beta,
                                                Double_t gamma,
                                                Double_t charge,
                                                Double_t p0,
                                                Double_t p1)
{
    // AZ-110124 - compute number of ionization clusters
 
    Double_t beta2 = beta * beta;
    Double_t gamma2 = gamma * gamma;
    Double_t val = p0 / beta2 * (p1 + TMath::Log(beta2 * gamma2) - beta2);
    return TMath::Min(val * charge * charge, 20000.);
}