diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc index a90bd123f7..fa27119df5 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc @@ -18,6 +18,8 @@ #include // for TrkrHitTruthA... #include #include +#include +#include #include #include @@ -52,26 +54,34 @@ #include #include +#include #include #include +#include #include +// for event id tagging in debug prints +#include + +#include + #include #include // for gsl_rng_alloc +#include #include #include #include // for sqrt, abs, NAN #include // for exit -#include #include #include // for _Rb_tree_cons... +#include #include // for pair namespace { template - constexpr T square(const T &x) + inline constexpr T square(const T &x) { return x * x; } @@ -108,8 +118,8 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl; exit(1); } - auto *runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); - auto *parNode = dynamic_cast(iter.findFirst("PHCompositeNode", "PAR")); + auto runNode = dynamic_cast(iter.findFirst("PHCompositeNode", "RUN")); + auto parNode = dynamic_cast(iter.findFirst("PHCompositeNode", "PAR")); const std::string paramnodename = "G4CELLPARAM_" + detector; const std::string geonodename = "G4CELLPAR_" + detector; const std::string tpcgeonodename = "G4GEO_" + detector; @@ -127,7 +137,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) if (!hitsetcontainer) { PHNodeIterator dstiter(dstNode); - auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); @@ -135,7 +145,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) } hitsetcontainer = new TrkrHitSetContainerv1; - auto *newNode = new PHIODataNode(hitsetcontainer, "TRKR_HITSET", "PHObject"); + auto newNode = new PHIODataNode(hitsetcontainer, "TRKR_HITSET", "PHObject"); DetNode->addNode(newNode); } @@ -143,7 +153,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) if (!hittruthassoc) { PHNodeIterator dstiter(dstNode); - auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); @@ -151,7 +161,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) } hittruthassoc = new TrkrHitTruthAssocv1; - auto *newNode = new PHIODataNode(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject"); + auto newNode = new PHIODataNode(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject"); DetNode->addNode(newNode); } @@ -159,7 +169,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) if (!truthtracks) { PHNodeIterator dstiter(dstNode); - auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); @@ -167,7 +177,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) } truthtracks = new TrkrTruthTrackContainerv1(); - auto *newNode = new PHIODataNode(truthtracks, "TRKR_TRUTHTRACKCONTAINER", "PHObject"); + auto newNode = new PHIODataNode(truthtracks, "TRKR_TRUTHTRACKCONTAINER", "PHObject"); DetNode->addNode(newNode); } @@ -175,7 +185,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) if (!truthclustercontainer) { PHNodeIterator dstiter(dstNode); - auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); @@ -183,22 +193,90 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) } truthclustercontainer = new TrkrClusterContainerv4; - auto *newNode = new PHIODataNode(truthclustercontainer, "TRKR_TRUTHCLUSTERCONTAINER", "PHObject"); + auto newNode = new PHIODataNode(truthclustercontainer, "TRKR_TRUTHCLUSTERCONTAINER", "PHObject"); DetNode->addNode(newNode); } - seggeonodename = "TPCGEOMCONTAINER"; // + detector; + m_track_map = findNode::getClass(topNode, "SvtxTrackMap"); + if (Verbosity() > 0 && m_track_map) + { + std::cout << Name() << ": found SvtxTrackMap with " << m_track_map->size() << " entries" << std::endl; + } + + seggeonodename = "CYLINDERCELLGEOM_SVTX"; // + detector; seggeo = findNode::getClass(topNode, seggeonodename); assert(seggeo); - // the z geometry is the same for all layers - PHG4TpcGeom *layergeom = seggeo->GetLayerCellGeom(20); - // from top of GEM stack to top of GEM stack - tpc_length = 2 * (layergeom->get_max_driftlength() + layergeom->get_CM_halfwidth()); - UpdateParametersWithMacro(); + + m_density_enabled = false; + m_density_layer = get_int_param("density_layer"); + m_density_bin_width_cm = get_double_param("density_bin_width_cm"); + m_density_window_cm = get_double_param("density_window_cm"); + m_uniform_density_test = (get_int_param("uniform_density_test") != 0); + if (m_density_layer >= 0 && m_density_bin_width_cm > 0.0 && m_density_window_cm > 0.0) + { + if (auto *layer_geom = seggeo->GetLayerCellGeom(m_density_layer)) + { + m_density_layer_radius = layer_geom->get_radius(); + const double half_thickness = 0.5 * layer_geom->get_thickness(); + m_density_layer_rlow = m_density_layer_radius - half_thickness; + m_density_layer_rhigh = m_density_layer_radius + half_thickness; + m_density_radial_margin_cm = std::max(0.0, get_double_param("density_radial_margin_cm")); + m_density_capture_rlow = m_density_layer_rlow - m_density_radial_margin_cm; + m_density_capture_rhigh = m_density_layer_rhigh + m_density_radial_margin_cm; + m_density_hist_half_range = (m_density_window_cm > 0.0) ? m_density_window_cm : (half_thickness + m_density_radial_margin_cm); + const double low_edge = -m_density_hist_half_range; + double high_edge = m_density_hist_half_range; + int nbins = static_cast(std::ceil((high_edge - low_edge) / m_density_bin_width_cm)); + if (nbins < 1) + { + nbins = 1; + } + high_edge = low_edge + nbins * m_density_bin_width_cm; + if (electronDensityProfile) + { + delete electronDensityProfile; + electronDensityProfile = nullptr; + } + const std::string hname = Name() + std::string("_ElectronDensityLayer") + std::to_string(m_density_layer); + const std::string htitle = "Ionization start vs. layer " + std::to_string(m_density_layer) + ";r - R_{layer} (cm);Electrons"; + electronDensityProfile = new TH1F(hname.c_str(), htitle.c_str(), nbins, low_edge, high_edge); + electronDensityProfile->SetDirectory(nullptr); + m_density_enabled = true; + } + else if (Verbosity() > 0) + { + std::cout << Name() << ": density_layer " << m_density_layer << " not present in geometry; disabling density histogram" << std::endl; + } + } + m_avg_x_enabled = false; + m_avg_layer = get_int_param("average_x_layer"); + if (m_avg_layer >= 0) + { + if (auto *layer_geom = seggeo->GetLayerCellGeom(m_avg_layer)) + { + m_avg_layer_radius = layer_geom->get_radius(); + const double half_thickness = 0.5 * layer_geom->get_thickness(); + m_avg_layer_rlow = m_avg_layer_radius - half_thickness; + m_avg_layer_rhigh = m_avg_layer_radius + half_thickness; + m_avg_x_enabled = true; + if (Verbosity() > 0) + { + const double full_thickness = layer_geom->get_thickness(); + std::cout << Name() << ": averaging x for layer " << m_avg_layer + << " at R=" << m_avg_layer_radius << " cm" + << " thickness=" << full_thickness << " cm (rlow=" << m_avg_layer_rlow + << ", rhigh=" << m_avg_layer_rhigh << ")" << std::endl; + } + } + else if (Verbosity() > 0) + { + std::cout << Name() << ": average_x_layer " << m_avg_layer << " not present in geometry; disabling average-x residual histogram" << std::endl; + } + } PHNodeIterator runIter(runNode); - auto *RunDetNode = dynamic_cast(runIter.findFirst("PHCompositeNode", detector)); + auto RunDetNode = dynamic_cast(runIter.findFirst("PHCompositeNode", detector)); if (!RunDetNode) { RunDetNode = new PHCompositeNode(detector); @@ -208,7 +286,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) // save this to the parNode for use PHNodeIterator parIter(parNode); - auto *ParDetNode = dynamic_cast(parIter.findFirst("PHCompositeNode", detector)); + auto ParDetNode = dynamic_cast(parIter.findFirst("PHCompositeNode", detector)); if (!ParDetNode) { ParDetNode = new PHCompositeNode(detector); @@ -218,11 +296,11 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) // find Tpc Geo PHNodeIterator tpcpariter(ParDetNode); - auto *tpcparams = findNode::getClass(ParDetNode, tpcgeonodename); + auto tpcparams = findNode::getClass(ParDetNode, tpcgeonodename); if (!tpcparams) { const std::string runparamname = "G4GEOPARAM_" + detector; - auto *tpcpdbparams = findNode::getClass(RunDetNode, runparamname); + auto tpcpdbparams = findNode::getClass(RunDetNode, runparamname); if (tpcpdbparams) { tpcparams = new PHParametersContainer(detector); @@ -247,6 +325,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) } const PHParameters *tpcparam = tpcparams->GetParameters(0); assert(tpcparam); + tpc_length = tpcparam->get_double_param("tpc_length"); diffusion_long = get_double_param("diffusion_long"); added_smear_sigma_long = get_double_param("added_smear_long"); @@ -256,30 +335,36 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) diffusion_trans *= zero_bfield_diffusion_factor; } added_smear_sigma_trans = get_double_param("added_smear_trans"); + drift_velocity = get_double_param("drift_velocity"); + force_min_trans_drift_length = get_double_param("force_min_trans_drift_length"); + if (force_min_trans_drift_length < 0.) + { + force_min_trans_drift_length = 0.; + } // Data on gasses @20 C and 760 Torr from the following source: // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations: // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html - double Ne_dEdx = 1.56; // keV/cm - double Ne_NTotal = 43; // Number/cm + double Ne_dEdx = 1.56; // keV/cm + double Ne_NTotal = 43; // Number/cm double Ne_frac = tpcparam->get_double_param("Ne_frac"); - double Ar_dEdx = 2.44; // keV/cm - double Ar_NTotal = 94; // Number/cm + double Ar_dEdx = 2.44; // keV/cm + double Ar_NTotal = 23; // Primary electrons/cm double Ar_frac = tpcparam->get_double_param("Ar_frac"); - double CF4_dEdx = 7; // keV/cm - double CF4_NTotal = 100; // Number/cm + double CF4_dEdx = 7; // keV/cm + double CF4_NTotal = 51; // Primary electrons/cm double CF4_frac = tpcparam->get_double_param("CF4_frac"); - double N2_dEdx = 2.127; // keV/cm https://pdg.lbl.gov/2024/AtomicNuclearProperties/HTML/nitrogen_gas.html - double N2_NTotal = 25; // Number/cm (probably not right but has a very small impact) + double N2_dEdx = 2.127; // keV/cm https://pdg.lbl.gov/2024/AtomicNuclearProperties/HTML/nitrogen_gas.html + double N2_NTotal = 25; // Number/cm (probably not right but has a very small impact) double N2_frac = tpcparam->get_double_param("N2_frac"); double isobutane_dEdx = 5.93; // keV/cm - double isobutane_NTotal = 195; // Number/cm + double isobutane_NTotal = 84; // Primary electrons/cm double isobutane_frac = tpcparam->get_double_param("isobutane_frac"); if (m_use_PDG_gas_params) @@ -291,39 +376,213 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) Ar_NTotal = 97; CF4_dEdx = 6.382; - CF4_NTotal = 120; + CF4_NTotal = 120; } - double Tpc_NTot = (Ne_NTotal * Ne_frac) + (Ar_NTotal * Ar_frac) + (CF4_NTotal * CF4_frac) + (N2_NTotal * N2_frac) + (isobutane_NTotal * isobutane_frac); + double Tpc_NTot = (Ne_NTotal * Ne_frac) + + (Ar_NTotal * Ar_frac) + + (CF4_NTotal * CF4_frac) + + (N2_NTotal * N2_frac) + + (isobutane_NTotal * isobutane_frac); + + double Tpc_dEdx = (Ne_dEdx * Ne_frac) + + (Ar_dEdx * Ar_frac) + + (CF4_dEdx * CF4_frac) + + (N2_dEdx * N2_frac) + + (isobutane_dEdx * isobutane_frac); + + electrons_per_gev = (Tpc_NTot / Tpc_dEdx) * 1e6; + + std::cout << "PHG4TpcElectronDrift::InitRun - electrons per GeV = " << electrons_per_gev << std::endl; + + // Initialize cluster size CDF for sPHENIX Ar/CF4/iC4H10 mixture. + constexpr int kClusterSizeCutoff = 384; + // Argon probabilities from the provided table (column b), converted to fractions. + const std::map ar_cluster_prob = { + {1, 0.656}, + {2, 0.150}, + {3, 0.064}, + {4, 0.035}, + {5, 0.0225}, + {6, 0.0155}, + {7, 0.0105}, + {8, 0.0081}, + {9, 0.0061}, + {10, 0.0049}, + {11, 0.0039}, + {12, 0.0030}, + {13, 0.0025}, + {14, 0.0020}, + {15, 0.0016}, + {16, 0.0012}, + {17, 0.00095}, + {18, 0.00075}, + {19, 0.00063}}; + // CH4 probabilities from the provided table (column b), converted to fractions. + const std::map ch4_cluster_prob = { + {1, 0.786}, + {2, 0.120}, + {3, 0.034}, + {4, 0.016}, + {5, 0.0095}, + {6, 0.0060}, + {7, 0.0044}, + {8, 0.0034}, + {9, 0.0027}, + {10, 0.0021}, + {11, 0.0017}, + {12, 0.0013}, + {13, 0.0010}, + {14, 0.0008}, + {15, 0.0006}, + {16, 0.00050}, + {17, 0.00042}, + {18, 0.00037}, + {19, 0.00033}}; + + auto build_cluster_prob = [&](const std::map &discrete, int tail_start) + { + std::vector prob(kClusterSizeCutoff + 1, 0.0); + double sum_discrete = 0.0; + for (auto const& [n, p] : discrete) + { + if (n > kClusterSizeCutoff) + { + continue; + } + prob[n] = p; + sum_discrete += p; + } + double remaining = 1.0 - sum_discrete; + if (remaining > 0.0 && tail_start <= kClusterSizeCutoff) + { + double sum_inv_sq = 0.0; + for (int n = tail_start; n <= kClusterSizeCutoff; ++n) + { + sum_inv_sq += 1.0 / (n * n); + } + double C = (sum_inv_sq > 0.0) ? remaining / sum_inv_sq : 0.0; + for (int n = tail_start; n <= kClusterSizeCutoff; ++n) + { + if (prob[n] == 0.0) + { + prob[n] = C / (n * n); + } + } + } + double norm = 0.0; + for (int n = 1; n <= kClusterSizeCutoff; ++n) + { + norm += prob[n]; + } + if (norm > 0.0) + { + for (int n = 1; n <= kClusterSizeCutoff; ++n) + { + prob[n] /= norm; + } + } + return prob; + }; + + const auto prob_ar = build_cluster_prob(ar_cluster_prob, 20); + const auto prob_ch4 = build_cluster_prob(ch4_cluster_prob, 20); + + const double mix_primary = (Ar_frac * Ar_NTotal) + + (CF4_frac * CF4_NTotal) + + (isobutane_frac * isobutane_NTotal); + const double weight_ar = (mix_primary > 0.0) ? (Ar_frac * Ar_NTotal / mix_primary) : 0.0; + const double weight_cf4 = (mix_primary > 0.0) ? (CF4_frac * CF4_NTotal / mix_primary) : 0.0; + const double weight_iso = (mix_primary > 0.0) ? (isobutane_frac * isobutane_NTotal / mix_primary) : 0.0; - double Tpc_dEdx = (Ne_dEdx * Ne_frac) + (Ar_dEdx * Ar_frac) + (CF4_dEdx * CF4_frac) + (N2_dEdx * N2_frac) + (isobutane_dEdx * isobutane_frac); + std::vector mix_prob(kClusterSizeCutoff + 1, 0.0); + if (mix_primary > 0.0) + { + for (int n = 1; n <= kClusterSizeCutoff; ++n) + { + mix_prob[n] = (weight_ar + weight_cf4) * prob_ar[n] + weight_iso * prob_ch4[n]; + } + } + else + { + mix_prob = prob_ar; + std::cout << "PHG4TpcElectronDrift::InitRun - Warning: gas mixture has no Ar, CF4, or isobutane content, defaulting to pure Argon cluster size distribution." << std::endl; + } - electrons_per_gev = (Tpc_NTot / Tpc_dEdx) * 1e6; + double mix_norm = 0.0; + for (int n = 1; n <= kClusterSizeCutoff; ++n) + { + mix_norm += mix_prob[n]; + } + if (mix_norm > 0.0) + { + for (int n = 1; n <= kClusterSizeCutoff; ++n) + { + mix_prob[n] /= mix_norm; + } + } - std::cout<<"PHG4ElectronDrift:: "<get_max_driftlength() / layergeom->get_drift_velocity_sim() + layergeom->get_extended_readout_time(); + max_time = get_double_param("max_time") + get_double_param("extended_readout_time"); min_active_radius = get_double_param("min_active_radius"); max_active_radius = get_double_param("max_active_radius"); if (Verbosity() > 0) { - std::cout << PHWHERE << " drift velocity " << layergeom->get_drift_velocity_sim() << " extended_readout_time " << layergeom->get_extended_readout_time() << " max time cutoff " << max_time << std::endl; + std::cout << Name() << " gas mixture fractions (Ne/Ar/CF4/N2/iC4H10): " + << Ne_frac << "/" << Ar_frac << "/" << CF4_frac << "/" << N2_frac << "/" << isobutane_frac << std::endl; + std::cout << Name() << " primary ionization summary: electrons per cm " << Tpc_NTot + << ", dE/dx (keV/cm) " << Tpc_dEdx + << ", electrons per GeV " << electrons_per_gev << std::endl; + std::cout << Name() << " diffusion sigmas (long/trans) [cm^0.5]: " << diffusion_long << "/" << diffusion_trans + << " with additional smearing (long/trans): " << added_smear_sigma_long << "/" << added_smear_sigma_trans << std::endl; + std::cout << Name() << " drift window [min,max] (ns): " << min_time << ", " << max_time << std::endl; + std::cout << Name() << " acceptance radii [min,max] (cm): " << min_active_radius << ", " << max_active_radius << std::endl; + if (force_min_trans_drift_length > 0.) + { + std::cout << Name() << " forcing minimum transverse drift length (cm): " << force_min_trans_drift_length << std::endl; + } + std::cout << PHWHERE << " drift velocity " << drift_velocity << " extended_readout_time " << get_double_param("extended_readout_time") << " max time cutoff " << max_time << std::endl; } - auto *se = Fun4AllServer::instance(); - dlong = new TH1F("difflong", "longitudinal diffusion", 100, diffusion_long - diffusion_long / 2., diffusion_long + diffusion_long / 2.); - se->registerHisto(dlong); - dtrans = new TH1F("difftrans", "transversal diffusion", 100, diffusion_trans - diffusion_trans / 2., diffusion_trans + diffusion_trans / 2.); - se->registerHisto(dtrans); + auto se = Fun4AllServer::instance(); + if (do_ElectronDriftQAHistos) + { + dlong = new TH1F("difflong", "longitudinal diffusion", 100, diffusion_long - diffusion_long / 2., diffusion_long + diffusion_long / 2.); + se->registerHisto(dlong); + dtrans = new TH1F("difftrans", "transversal diffusion", 100, diffusion_trans - diffusion_trans / 2., diffusion_trans + diffusion_trans / 2.); + se->registerHisto(dtrans); + } - do_ElectronDriftQAHistos = false; // Whether or not to produce an ElectronDriftQA.root file with useful info if (do_ElectronDriftQAHistos) { + diffDistance = new TH1F("diffDistance", "Transverse diffusion displacement;#Delta r (cm);Counts", 300, 0.0, 3.0); + diffDX = new TH1F("diffDX", "Transverse diffusion #Delta x;#Delta x (cm);Counts", 400, -3.0, 3.0); + diffDY = new TH1F("diffDY", "Transverse diffusion #Delta y;#Delta y (cm);Counts", 400, -3.0, 3.0); + nElectrons = new TH1F("nElectrons", "Sampled electrons per step;N_{e};Counts", 400, -0.5, 399.5); + poissonMean = new TH1F("electronMean", "Poisson mean per step;#bar{N}_{e};Counts", 400, 0.0, 200.0); + nElectronsVsMean = new TH2F("nElectronsVsMean", "Sampled electrons vs. mean;#bar{N}_{e};N_{e}", 200, 0.0, 200.0, 200, -0.5, 199.5); + diffPerSqrtL = new TH1F("diffPerSqrtL", "Transverse diffusion normalized by #sqrt{L};#Delta r/#sqrt{L} (cm^{0.5});Counts", 400, 0.0, 0.2); + diffDXPerSqrtL = new TH1F("diffDXPerSqrtL", "#Delta x normalized by #sqrt{L};#Delta x/#sqrt{L} (cm^{0.5});Counts", 400, -0.2, 0.2); + diffDYPerSqrtL = new TH1F("diffDYPerSqrtL", "#Delta y normalized by #sqrt{L};#Delta y/#sqrt{L} (cm^{0.5});Counts", 400, -0.2, 0.2); + nElectronsPerCm = new TH1F("nElectronsPerCm", "Sampled electrons per cm;N_{e}/cm;Counts", 400, 0.0, 400.0); + diffVsDrift = new TH2F("diffVsDrift", "Transverse diffusion vs. drift length;Drift length L (cm);#Delta r (cm)", 200, 0.0, tpc_length / 2., 300, 0.0, 3.0); + diffDXVsDrift = new TH2F("diffDXVsDrift", "#Delta x vs. drift length;Drift length L (cm);#Delta x (cm)", 200, 0.0, tpc_length / 2., 400, -3.0, 3.0); + diffDYVsDrift = new TH2F("diffDYVsDrift", "#Delta y vs. drift length;Drift length L (cm);#Delta y (cm)", 200, 0.0, tpc_length / 2., 400, -3.0, 3.0); + diffPerSqrtLVsDrift = new TH2F("diffPerSqrtLVsDrift", "#Delta r/#sqrt{L} vs. drift length;Drift length L (cm);#Delta r/#sqrt{L} (cm^{0.5})", 200, 0.0, tpc_length / 2., 300, 0.0, 0.2); + diffDXPerSqrtLVsDrift = new TH2F("diffDXPerSqrtLVsDrift", "#Delta x/#sqrt{L} vs. drift length;Drift length L (cm);#Delta x/#sqrt{L} (cm^{0.5})", 200, 0.0, tpc_length / 2., 400, -0.2, 0.2); + diffDYPerSqrtLVsDrift = new TH2F("diffDYPerSqrtLVsDrift", "#Delta y/#sqrt{L} vs. drift length;Drift length L (cm);#Delta y/#sqrt{L} (cm^{0.5})", 200, 0.0, tpc_length / 2., 400, -0.2, 0.2); hitmapstart = new TH2F("hitmapstart", "g4hit starting X-Y locations", 1560, -78, 78, 1560, -78, 78); hitmapend = new TH2F("hitmapend", "g4hit final X-Y locations", 1560, -78, 78, 1560, -78, 78); hitmapstart_z = new TH2F("hitmapstart_z", "g4hit starting Z-R locations", 2000, -100, 100, 780, 0, 78); @@ -339,6 +598,42 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) deltarnodiff = new TH2F("deltarnodiff", "Delta r (no diffusion, only SC distortion); r (cm);#Delta r (cm)", 580, 20, 78, 1000, -2, 5); deltarnodist = new TH2F("deltarnodist", "Delta r (no SC distortion, only diffusion); r (cm);#Delta r (cm)", 580, 20, 78, 1000, -2, 5); ratioElectronsRR = new TH1F("ratioElectronsRR", "Ratio of electrons reach readout vs all in acceptance", 1561, -0.0325, 1.0465); + driftXY = new TNtuple("driftXY", "Electron start/end positions", "xs:ys:xf:yf"); + } + + if (m_avg_x_enabled) + { + m_avgOutf.reset(new TFile(m_avg_output_file.c_str(), "recreate")); + if (!m_avgOutf || m_avgOutf->IsZombie()) + { + std::cout << Name() << ": failed to open avg output file " << m_avg_output_file << std::endl; + } + if (m_avgOutf) + { + m_avgOutf->cd(); + } + m_avgXResidualTree = new TTree("avgXResidual", "Track-level -x_int residuals"); + m_avgXResidualTree->Branch("trackid", &m_avgTree_trackid, "trackid/F"); + m_avgXResidualTree->Branch("residual", &m_avgTree_residual, "residual/F"); + m_avgXResidualTree->Branch("avgx", &m_avgTree_avgx, "avgx/F"); + m_avgXResidualTree->Branch("xint", &m_avgTree_xint, "xint/F"); + m_avgXResidualTree->Branch("residual_primary", &m_avgTree_residual_primary, "residual_primary/F"); + m_avgXResidualTree->Branch("avgx_primary", &m_avgTree_avgx_primary, "avgx_primary/F"); + m_avgXResidualTree->Branch("count_primary", &m_avgTree_count_primary, "count_primary/F"); + m_avgXResidualTree->Branch("residual_end", &m_avgTree_residual_end, "residual_end/F"); + m_avgXResidualTree->Branch("avgx_end", &m_avgTree_avgx_end, "avgx_end/F"); + m_avgXResidualTree->Branch("count_end", &m_avgTree_count_end, "count_end/F"); + m_avgXResidualTree->Branch("residual_end_primary", &m_avgTree_residual_end_primary, "residual_end_primary/F"); + m_avgXResidualTree->Branch("avgx_end_primary", &m_avgTree_avgx_end_primary, "avgx_end_primary/F"); + m_avgXResidualTree->Branch("count_end_primary", &m_avgTree_count_end_primary, "count_end_primary/F"); + m_avgXResidualTree->Branch("residual_x", &m_avgTree_residual_x, "residual_x/F"); + m_avgXResidualTree->Branch("avgx_cart", &m_avgTree_avgx_cart, "avgx_cart/F"); + m_avgXResidualTree->Branch("residual_x_primary", &m_avgTree_residual_x_primary, "residual_x_primary/F"); + m_avgXResidualTree->Branch("avgx_cart_primary", &m_avgTree_avgx_cart_primary, "avgx_cart_primary/F"); + m_avgXResidualTree->Branch("residual_x_end", &m_avgTree_residual_x_end, "residual_x_end/F"); + m_avgXResidualTree->Branch("avgx_cart_end", &m_avgTree_avgx_cart_end, "avgx_cart_end/F"); + m_avgXResidualTree->Branch("residual_x_end_primary", &m_avgTree_residual_x_end_primary, "residual_x_end_primary/F"); + m_avgXResidualTree->Branch("avgx_cart_end_primary", &m_avgTree_avgx_cart_end_primary, "avgx_cart_end_primary/F"); } if (Verbosity()) @@ -365,7 +660,7 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) for (auto layeriter = range.first; layeriter != range.second; ++layeriter) { const auto radius = layeriter->second->get_radius(); - std::cout << std::format("{:.3f} ", radius); + std::cout << boost::str(boost::format("%.3f ") % radius); if (++counter == 8) { counter = 0; @@ -382,14 +677,14 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) if (!mClusHitsVerbose) { PHNodeIterator dstiter(dstNode); - auto *DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); + auto DetNode = dynamic_cast(dstiter.findFirst("PHCompositeNode", "TRKR")); if (!DetNode) { DetNode = new PHCompositeNode("TRKR"); dstNode->addNode(DetNode); } mClusHitsVerbose = new ClusHitsVerbosev1(); - auto *newNode = new PHIODataNode(mClusHitsVerbose, "Trkr_TruthClusHitsVerbose", "PHObject"); + auto newNode = new PHIODataNode(mClusHitsVerbose, "Trkr_TruthClusHitsVerbose", "PHObject"); DetNode->addNode(newNode); } } @@ -407,16 +702,15 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) std::cout << PHWHERE << "ActsGeometry not found on node tree. Exiting" << std::endl; return Fun4AllReturnCodes::ABORTRUN; } + m_track_layer_data.clear(); - PHG4TpcGeom *layergeom = seggeo->GetLayerCellGeom(20); - if (truth_clusterer.needs_input_nodes()) { truth_clusterer.set_input_nodes(truthclustercontainer, m_tGeometry, seggeo, mClusHitsVerbose); } - static constexpr unsigned int print_layer = 18; + static constexpr unsigned int print_layer = 49; // tells m_distortionMap which event to look at if (m_distortionMap) @@ -425,7 +719,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) } // g4hits - auto *g4hit = findNode::getClass(topNode, hitnodename); + auto g4hit = findNode::getClass(topNode, hitnodename); if (!g4hit) { std::cout << "Could not locate g4 hit node " << hitnodename << std::endl; @@ -438,8 +732,12 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) unsigned int count_g4hits = 0; // int count_electrons = 0; + m_track_path_offset.clear(); + m_track_anchor_set.clear(); + m_track_anchor_length.clear(); + // double ecollectedhits = 0.0; - // int ncollectedhits = 0; +// int ncollectedhits = 0; double ihit = 0; unsigned int dump_interval = 5000; // dump temp_hitsetcontainer to the node tree after this many g4hits unsigned int dump_counter = 0; @@ -494,6 +792,30 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) std::cout << " NOT embedded" << std::endl; } } + + if (m_avg_x_enabled && m_track_map && trkid >= 0) + { + auto it = m_track_map->find(static_cast(trkid)); + if (it != m_track_map->end()) + { + const SvtxTrack *trk = it->second; + const double vx = trk->get_px(); + const double vy = trk->get_py(); + const double vz = trk->get_pz(); + const double vmag = std::sqrt(square(vx) + square(vy) + square(vz)); + if (vmag > 0.) + { + auto &layer_data = m_track_layer_data[trkid]; + layer_data.base_x = trk->get_x(); + layer_data.base_y = trk->get_y(); + layer_data.base_z = trk->get_z(); + layer_data.dir_x = vx / vmag; + layer_data.dir_y = vy / vmag; + layer_data.dir_z = vz / vmag; + layer_data.have_line = true; + } + } + } } // see if there is a jump in x or y relative to previous PHG4Hit @@ -519,23 +841,91 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) // drifted electrons, then copy to the node tree later double eion = hiter->second->get_eion(); - unsigned int n_electrons = gsl_ran_poisson(RandomGenerator.get(), eion * electrons_per_gev); + const double poisson_mean = eion * electrons_per_gev; + unsigned int n_electrons = 0; + if (m_uniform_density_test) + { + n_electrons = (poisson_mean > 0.) ? static_cast(std::lround(poisson_mean)) : 0; + } + else + { + n_electrons = gsl_ran_poisson(RandomGenerator.get(), poisson_mean); + } // count_electrons += n_electrons; + const double dx_hit = hiter->second->get_x(1) - hiter->second->get_x(0); + const double dy_hit = hiter->second->get_y(1) - hiter->second->get_y(0); + const double dz_hit = hiter->second->get_z(1) - hiter->second->get_z(0); + const double step_length = std::sqrt(square(dx_hit) + square(dy_hit) + square(dz_hit)); + if (m_avg_x_enabled) + { + auto &layer_data = m_track_layer_data[trkid_new]; + if (!layer_data.have_line) + { + const double norm = step_length; + if (norm > 0.) + { + layer_data.base_x = hiter->second->get_x(0); + layer_data.base_y = hiter->second->get_y(0); + layer_data.base_z = hiter->second->get_z(0); + layer_data.dir_x = dx_hit / norm; + layer_data.dir_y = dy_hit / norm; + layer_data.dir_z = dz_hit / norm; + layer_data.have_line = true; + } + } + } + + if (do_ElectronDriftQAHistos) + { + if (poissonMean) + { + poissonMean->Fill(poisson_mean); + } + if (nElectrons) + { + nElectrons->Fill(static_cast(n_electrons)); + } + if (nElectronsPerCm && step_length > 0.) + { + nElectronsPerCm->Fill(static_cast(n_electrons) / step_length); + } + if (nElectronsVsMean) + { + nElectronsVsMean->Fill(poisson_mean, static_cast(n_electrons)); + } + } + if (Verbosity() > 100) { std::cout << " new hit with t0, " << t0 << " g4hitid " << hiter->first - << " eion " << eion << " n_electrons " << n_electrons + << " eion " << eion << " poisson mean " << poisson_mean + << " n_electrons " << n_electrons << " entry z " << hiter->second->get_z(0) << " exit z " << hiter->second->get_z(1) << " avg z" << (hiter->second->get_z(0) + hiter->second->get_z(1)) / 2.0 << std::endl; } + if (Verbosity() > 0 && count_g4hits <= 5) + { + std::cout << Name() << " hit " << count_g4hits + << " eion " << eion + << " poisson mean " << poisson_mean + << " sampled electrons " << n_electrons + << " track dz " << hiter->second->get_z(1) - hiter->second->get_z(0) + << " entry radius " << std::sqrt(square(hiter->second->get_x(0)) + square(hiter->second->get_y(0))) + << " exit radius " << std::sqrt(square(hiter->second->get_x(1)) + square(hiter->second->get_y(1))) + << std::endl; + } + if (n_electrons == 0) { + m_track_path_offset[hiter->second->get_trkid()] += step_length; continue; } + const int track_id = hiter->second->get_trkid(); + double track_segment_start = m_track_path_offset[track_id]; if (Verbosity() > 100) { @@ -551,36 +941,68 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) } int notReachingReadout = 0; - // int notInAcceptance = 0; +// int notInAcceptance = 0; + + // Loop over primary electrons (clusters) for (unsigned int i = 0; i < n_electrons; i++) { + // Sample cluster size + // Sample cluster size + int cluster_size = 1; + if (m_enable_laser_clustering) + { + double p = gsl_ran_flat(RandomGenerator.get(), 0.0, 1.0); + auto it = std::lower_bound(cluster_size_cdf.begin(), cluster_size_cdf.end(), p); + if (it != cluster_size_cdf.end()) + { + cluster_size = std::distance(cluster_size_cdf.begin(), it); + } + } + // We choose the electron starting position at random from a flat // distribution along the path length the parameter t is the fraction of // the distance along the path betwen entry and exit points, it has // values between 0 and 1 - const double f = gsl_ran_flat(RandomGenerator.get(), 0.0, 1.0); + double f = gsl_ran_flat(RandomGenerator.get(), 0.0, 1.0); + if (m_uniform_density_test && n_electrons > 0) + { + f = (static_cast(i) + 0.5) / static_cast(n_electrons); + if (f >= 1.0) + { + f = std::nextafter(1.0, 0.0); + } + } + + // Loop over secondary electrons in the cluster + for (int j = 0; j < cluster_size; ++j) + { + const double x_start = hiter->second->get_x(0) + f * (hiter->second->get_x(1) - hiter->second->get_x(0)); const double y_start = hiter->second->get_y(0) + f * (hiter->second->get_y(1) - hiter->second->get_y(0)); const double z_start = hiter->second->get_z(0) + f * (hiter->second->get_z(1) - hiter->second->get_z(0)); const double t_start = hiter->second->get_t(0) + f * (hiter->second->get_t(1) - hiter->second->get_t(0)); - unsigned int side = 0; if (z_start > 0) { side = 1; } - const double r_sigma = diffusion_trans * sqrt(tpc_length / 2. - std::abs(z_start)); - const double rantrans = + double drift_distance = tpc_length / 2. - std::abs(z_start); + if (drift_distance < 0.) + { + drift_distance = 0.; + } + const double r_sigma = diffusion_trans * sqrt(drift_distance); + double rantrans = gsl_ran_gaussian(RandomGenerator.get(), r_sigma) + gsl_ran_gaussian(RandomGenerator.get(), added_smear_sigma_trans); - const double t_path = (tpc_length / 2. - std::abs(z_start)) / layergeom->get_drift_velocity_sim(); - const double t_sigma = diffusion_long * sqrt(tpc_length / 2. - std::abs(z_start)) / layergeom->get_drift_velocity_sim(); + const double t_path = drift_distance / drift_velocity; + const double t_sigma = diffusion_long * std::sqrt(drift_distance) / drift_velocity; const double rantime = gsl_ran_gaussian(RandomGenerator.get(), t_sigma) + - gsl_ran_gaussian(RandomGenerator.get(), added_smear_sigma_long) / layergeom->get_drift_velocity_sim(); + gsl_ran_gaussian(RandomGenerator.get(), added_smear_sigma_long) / drift_velocity; double t_final = t_start + t_path + rantime; if (t_final < min_time || t_final > max_time) @@ -591,25 +1013,126 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) double z_final; if (z_start < 0) { - z_final = -tpc_length / 2. + t_final * layergeom->get_drift_velocity_sim(); + z_final = -tpc_length / 2. + t_final * drift_velocity; } else { - z_final = tpc_length / 2. - t_final * layergeom->get_drift_velocity_sim(); + z_final = tpc_length / 2. - t_final * drift_velocity; } const double radstart = std::sqrt(square(x_start) + square(y_start)); const double phistart = std::atan2(y_start, x_start); + if (m_avg_x_enabled && radstart >= m_avg_layer_rlow && radstart <= m_avg_layer_rhigh) + { + // collect undiffused electron start positions for the configured layer + auto &layer_data = m_track_layer_data[track_id]; + layer_data.sum_rphi += m_avg_layer_radius * phistart; + layer_data.sum_x += x_start; + layer_data.count++; + if (j == 0) + { + layer_data.sum_rphi_primary += m_avg_layer_radius * phistart; + layer_data.sum_x_primary += x_start; + layer_data.count_primary++; + } + } const double ranphi = gsl_ran_flat(RandomGenerator.get(), -M_PI, M_PI); - double x_final = x_start + rantrans * std::cos(ranphi); // Initialize these to be only diffused first, will be overwritten if doing SC distortion double y_final = y_start + rantrans * std::sin(ranphi); + if (force_min_trans_drift_length > drift_distance) + { + const double additional_length = force_min_trans_drift_length - drift_distance; + if (additional_length > 0.) + { + const double extra_sigma = diffusion_trans * std::sqrt(additional_length); + const double extra_r = gsl_ran_gaussian(RandomGenerator.get(), extra_sigma); + const double extra_phi = gsl_ran_flat(RandomGenerator.get(), -M_PI, M_PI); + x_final += extra_r * std::cos(extra_phi); + y_final += extra_r * std::sin(extra_phi); + } + } + + const double delta_x = x_final - x_start; + const double delta_y = y_final - y_start; + rantrans = std::sqrt(square(delta_x) + square(delta_y)); + double rad_final = sqrt(square(x_final) + square(y_final)); double phi_final = atan2(y_final, x_final); + if (m_avg_x_enabled && rad_final >= m_avg_layer_rlow && rad_final <= m_avg_layer_rhigh) + { + auto &layer_data = m_track_layer_data[track_id]; + // Use the nominal layer radius so residuals are purely angular and not + // inflated by radius fluctuations from diffusion. + layer_data.sum_rphi_end += m_avg_layer_radius * phi_final; + layer_data.sum_x_end += x_final; + layer_data.count_end++; + if (j == 0) + { + layer_data.sum_rphi_end_primary += m_avg_layer_radius * phi_final; + layer_data.sum_x_end_primary += x_final; + layer_data.count_end_primary++; + } + } if (do_ElectronDriftQAHistos) { + if (diffDistance) + { + diffDistance->Fill(rantrans); + } + if (diffDX) + { + diffDX->Fill(delta_x); + } + if (diffDY) + { + diffDY->Fill(delta_y); + } + if (diffVsDrift) + { + diffVsDrift->Fill(drift_distance, rantrans); + } + if (diffDXVsDrift) + { + diffDXVsDrift->Fill(drift_distance, delta_x); + } + if (diffDYVsDrift) + { + diffDYVsDrift->Fill(drift_distance, delta_y); + } + + if (drift_distance > 0.) + { + const double sqrtL = std::sqrt(drift_distance); + const double norm = rantrans / sqrtL; + if (diffPerSqrtL) + { + diffPerSqrtL->Fill(norm); + } + if (diffPerSqrtLVsDrift) + { + diffPerSqrtLVsDrift->Fill(drift_distance, norm); + } + const double dxnorm = delta_x / sqrtL; + const double dynorm = delta_y / sqrtL; + if (diffDXPerSqrtL) + { + diffDXPerSqrtL->Fill(dxnorm); + } + if (diffDYPerSqrtL) + { + diffDYPerSqrtL->Fill(dynorm); + } + if (diffDXPerSqrtLVsDrift) + { + diffDXPerSqrtLVsDrift->Fill(drift_distance, dxnorm); + } + if (diffDYPerSqrtLVsDrift) + { + diffDYPerSqrtLVsDrift->Fill(drift_distance, dynorm); + } + } z_startmap->Fill(z_start, radstart); // map of starting location in Z vs. R deltaphinodist->Fill(phistart, rantrans / rad_final); // delta phi no distortion, just diffusion+smear deltarnodist->Fill(radstart, rantrans); // delta r no distortion, just diffusion+smear @@ -634,11 +1157,11 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) z_final += z_distortion; if (z_start < 0) { - t_final = (z_final + tpc_length / 2.0) / layergeom->get_drift_velocity_sim(); + t_final = (z_final + tpc_length / 2.0) / drift_velocity; } else { - t_final = (tpc_length / 2.0 - z_final) / layergeom->get_drift_velocity_sim(); + t_final = (tpc_length / 2.0 - z_final) / drift_velocity; } x_final = rad_final * std::cos(phi_final); @@ -667,10 +1190,33 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) } } + if (do_ElectronDriftQAHistos && driftXY) + { + driftXY->Fill(x_start, y_start, x_final, y_final); + } + + if (m_density_enabled && electronDensityProfile) + { + if (rad_final >= m_density_capture_rlow && rad_final <= m_density_capture_rhigh) + { + const double delta_r = radstart - m_density_layer_radius; + if (delta_r >= -m_density_hist_half_range && delta_r <= m_density_hist_half_range) + { + const double s_global = track_segment_start + f * step_length; + if (!m_track_anchor_set[track_id]) + { + m_track_anchor_set[track_id] = true; + m_track_anchor_length[track_id] = s_global; + } + electronDensityProfile->Fill(s_global - m_track_anchor_length[track_id]); + } + } + } + // remove electrons outside of our acceptance. Careful though, electrons from just inside 30 cm can contribute in the 1st active layer readout, so leave a little margin if (rad_final < min_active_radius - 2.0 || rad_final > max_active_radius + 1.0) { - // notInAcceptance++; +// notInAcceptance++; continue; } @@ -701,8 +1247,11 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) padplane->MapToPadPlane(truth_clusterer, single_hitsetcontainer.get(), temp_hitsetcontainer.get(), hittruthassoc, x_final, y_final, t_final, side, hiter, ntpad, nthit); + } // end loop over secondary electrons } // end loop over electrons for this g4hit + m_track_path_offset[track_id] = track_segment_start + step_length; + if (do_ElectronDriftQAHistos) { ratioElectronsRR->Fill((double) (n_electrons - notReachingReadout) / n_electrons); @@ -746,7 +1295,27 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) // - if this is the last g4hit if (dump_counter >= dump_interval || count_g4hits == g4hit->size()) { - // std::cout << " dump_counter " << dump_counter << " count_g4hits " << count_g4hits << std::endl; + // additional debug: entering dump/copy of temp_hitsetcontainer to node TRKR_HITSET + if (Verbosity() > 50) + { + int evtseq_dbg = -1; + if (auto* eh = findNode::getClass(topNode, "EventHeader")) + { + evtseq_dbg = eh->get_EvtSequence(); + } + // count how many temp hitsets we are about to copy + unsigned int ntemp = 0; + { + auto tmp_rng_dbg = temp_hitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId); + for (auto itdbg = tmp_rng_dbg.first; itdbg != tmp_rng_dbg.second; ++itdbg) ++ntemp; + } + std::cout << "EDrift: entering dump: evt=" << evtseq_dbg + << " dump_counter=" << dump_counter + << " count_g4hits=" << count_g4hits + << " g4hit_size=" << g4hit->size() + << " temp_hitsets=" << ntemp + << std::endl; + } double eg4hit = 0.0; TrkrHitSetContainer::ConstRange temp_hitset_range = temp_hitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId); @@ -759,7 +1328,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) const unsigned int layer = TrkrDefs::getLayer(node_hitsetkey); const int sector = TpcDefs::getSectorId(node_hitsetkey); const int side = TpcDefs::getSide(node_hitsetkey); - if (Verbosity() > 100) + if (Verbosity() > 50) { std::cout << "PHG4TpcElectronDrift: temp_hitset with key: " << node_hitsetkey << " in layer " << layer << " with sector " << sector << " side " << side << std::endl; @@ -784,7 +1353,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) eg4hit += temp_tpchit->getEnergy(); // ecollectedhits += temp_tpchit->getEnergy(); - // ncollectedhits++; +// ncollectedhits++; } // find or add this hit to the node tree @@ -801,7 +1370,7 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) } // end loop over temp hits - if (Verbosity() > 100 && layer == print_layer) + if (Verbosity() > 50 && layer == print_layer) { std::cout << " ihit " << ihit << " collected energy = " << eg4hit << std::endl; } @@ -821,6 +1390,146 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) } // end loop over g4hits + // Force dump any remaining hits that didn't get dumped in the loop + TrkrHitSetContainer::ConstRange temp_hitset_range_final = temp_hitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId); + if (temp_hitset_range_final.first != temp_hitset_range_final.second) // if there are any temp hitsets + { + if (Verbosity() > 50) + { + std::cout << "PHG4TpcElectronDrift: Forcing final dump of temp_hitsetcontainer" << std::endl; + } + + for (TrkrHitSetContainer::ConstIterator temp_hitset_iter = temp_hitset_range_final.first; + temp_hitset_iter != temp_hitset_range_final.second; + ++temp_hitset_iter) + { + TrkrDefs::hitsetkey node_hitsetkey = temp_hitset_iter->first; + const unsigned int layer = TrkrDefs::getLayer(node_hitsetkey); + const int sector = TpcDefs::getSectorId(node_hitsetkey); + const int side = TpcDefs::getSide(node_hitsetkey); + + if (Verbosity() > 50) + { + std::cout << "PHG4TpcElectronDrift: final dump - temp_hitset with key: " << node_hitsetkey + << " in layer " << layer << " with sector " << sector << " side " << side << std::endl; + } + + // find or add this hitset on the node tree + TrkrHitSetContainer::Iterator node_hitsetit = hitsetcontainer->findOrAddHitSet(node_hitsetkey); + + // get all of the hits from the temporary hitset + TrkrHitSet::ConstRange temp_hit_range = temp_hitset_iter->second->getHits(); + for (TrkrHitSet::ConstIterator temp_hit_iter = temp_hit_range.first; + temp_hit_iter != temp_hit_range.second; + ++temp_hit_iter) + { + TrkrDefs::hitkey temp_hitkey = temp_hit_iter->first; + TrkrHit *temp_tpchit = temp_hit_iter->second; + + // find or add this hit to the node tree + TrkrHit *node_hit = node_hitsetit->second->getHit(temp_hitkey); + if (!node_hit) + { + // Otherwise, create a new one + node_hit = new TrkrHitv2(); + node_hitsetit->second->addHitSpecificKey(temp_hitkey, node_hit); + } + + // Either way, add the energy to it + node_hit->addEnergy(temp_tpchit->getEnergy()); + } // end loop over temp hits + } // end loop over temp hitsets + + // erase all entries in the temp hitsetcontainer + temp_hitsetcontainer->Reset(); + } + + if (m_avg_x_enabled && m_avgXResidualTree) + { + for (const auto &entry : m_track_layer_data) + { + const auto &data = entry.second; + if (!data.have_line || data.count == 0) + { + continue; + } + double t_out = 0.0; + double xi = 0.0; + double yi = 0.0; + double zi = 0.0; + if (!cylinder_intersection(data.base_x, data.base_y, data.base_z, + data.dir_x, data.dir_y, data.dir_z, + m_avg_layer_radius, t_out, xi, yi, zi)) + { + continue; + } + const double avg_rphi = data.sum_rphi / static_cast(data.count); + const double rphi_int = m_avg_layer_radius * std::atan2(yi, xi); + m_avgTree_trackid = static_cast(entry.first); + m_avgTree_residual = static_cast(avg_rphi - rphi_int); + m_avgTree_avgx = static_cast(avg_rphi); + m_avgTree_xint = static_cast(rphi_int); + const double avg_x_cart = data.sum_x / static_cast(data.count); + m_avgTree_residual_x = static_cast(avg_x_cart - xi); + m_avgTree_avgx_cart = static_cast(avg_x_cart); + if (data.count_primary > 0) + { + const double avg_rphi_primary = data.sum_rphi_primary / static_cast(data.count_primary); + m_avgTree_residual_primary = static_cast(avg_rphi_primary - rphi_int); + m_avgTree_avgx_primary = static_cast(avg_rphi_primary); + m_avgTree_count_primary = static_cast(data.count_primary); + const double avg_x_primary = data.sum_x_primary / static_cast(data.count_primary); + m_avgTree_residual_x_primary = static_cast(avg_x_primary - xi); + m_avgTree_avgx_cart_primary = static_cast(avg_x_primary); + } + else + { + m_avgTree_residual_primary = std::numeric_limits::quiet_NaN(); + m_avgTree_avgx_primary = std::numeric_limits::quiet_NaN(); + m_avgTree_count_primary = 0.0f; + m_avgTree_residual_x_primary = std::numeric_limits::quiet_NaN(); + m_avgTree_avgx_cart_primary = std::numeric_limits::quiet_NaN(); + } + if (data.count_end > 0) + { + const double avg_rphi_end = data.sum_rphi_end / static_cast(data.count_end); + m_avgTree_residual_end = static_cast(avg_rphi_end - rphi_int); + m_avgTree_avgx_end = static_cast(avg_rphi_end); + m_avgTree_count_end = static_cast(data.count_end); + const double avg_x_end = data.sum_x_end / static_cast(data.count_end); + m_avgTree_residual_x_end = static_cast(avg_x_end - xi); + m_avgTree_avgx_cart_end = static_cast(avg_x_end); + } + else + { + m_avgTree_residual_end = std::numeric_limits::quiet_NaN(); + m_avgTree_avgx_end = std::numeric_limits::quiet_NaN(); + m_avgTree_count_end = 0.0f; + m_avgTree_residual_x_end = std::numeric_limits::quiet_NaN(); + m_avgTree_avgx_cart_end = std::numeric_limits::quiet_NaN(); + } + if (data.count_end_primary > 0) + { + const double avg_rphi_end_primary = data.sum_rphi_end_primary / static_cast(data.count_end_primary); + m_avgTree_residual_end_primary = static_cast(avg_rphi_end_primary - rphi_int); + m_avgTree_avgx_end_primary = static_cast(avg_rphi_end_primary); + m_avgTree_count_end_primary = static_cast(data.count_end_primary); + const double avg_x_end_primary = data.sum_x_end_primary / static_cast(data.count_end_primary); + m_avgTree_residual_x_end_primary = static_cast(avg_x_end_primary - xi); + m_avgTree_avgx_cart_end_primary = static_cast(avg_x_end_primary); + } + else + { + m_avgTree_residual_end_primary = std::numeric_limits::quiet_NaN(); + m_avgTree_avgx_end_primary = std::numeric_limits::quiet_NaN(); + m_avgTree_count_end_primary = 0.0f; + m_avgTree_residual_x_end_primary = std::numeric_limits::quiet_NaN(); + m_avgTree_avgx_cart_end_primary = std::numeric_limits::quiet_NaN(); + } + m_avgXResidualTree->Fill(); + } + } + if (truth_track) { truth_clusterer.cluster_hits(truth_track); @@ -829,7 +1538,9 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) if (Verbosity() > 20) { - std::cout << "From PHG4TpcElectronDrift: hitsetcontainer printout at end:" << std::endl; + int evtseq = -1; + if (auto* eh = findNode::getClass(topNode, "EventHeader")) { evtseq = eh->get_EvtSequence(); } + std::cout << "From PHG4TpcElectronDrift: hitsetcontainer printout at end: evt=" << evtseq << std::endl; // We want all hitsets for the Tpc TrkrHitSetContainer::ConstRange hitset_range = hitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId); for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first; @@ -839,10 +1550,10 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) // we have an itrator to one TrkrHitSet for the Tpc from the trkrHitSetContainer TrkrDefs::hitsetkey hitsetkey = hitset_iter->first; const unsigned int layer = TrkrDefs::getLayer(hitsetkey); - if (layer != print_layer) - { - continue; - } + //if (layer != print_layer) + //{ + // continue; + // } const int sector = TpcDefs::getSectorId(hitsetkey); const int side = TpcDefs::getSide(hitsetkey); @@ -851,6 +1562,8 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) // get all of the hits from this hitset TrkrHitSet *hitset = hitset_iter->second; TrkrHitSet::ConstRange hit_range = hitset->getHits(); + + for (TrkrHitSet::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter) @@ -907,7 +1620,7 @@ int PHG4TpcElectronDrift::End(PHCompositeNode * /*topNode*/) } if (do_ElectronDriftQAHistos) { - EDrift_outf.reset(new TFile("ElectronDriftQA.root", "recreate")); + EDrift_outf.reset(new TFile(m_qa_output_file.c_str(), "recreate")); EDrift_outf->cd(); deltar->Write(); deltaphi->Write(); @@ -924,8 +1637,90 @@ int PHG4TpcElectronDrift::End(PHCompositeNode * /*topNode*/) hitmapend_z->Write(); z_startmap->Write(); ratioElectronsRR->Write(); + if (diffDistance) + { + diffDistance->Write(); + } + if (diffDX) + { + diffDX->Write(); + } + if (diffDY) + { + diffDY->Write(); + } + if (diffPerSqrtL) + { + diffPerSqrtL->Write(); + } + if (nElectronsPerCm) + { + nElectronsPerCm->Write(); + } + if (diffDXPerSqrtL) + { + diffDXPerSqrtL->Write(); + } + if (diffDYPerSqrtL) + { + diffDYPerSqrtL->Write(); + } + if (nElectrons) + { + nElectrons->Write(); + } + if (poissonMean) + { + poissonMean->Write(); + } + if (nElectronsVsMean) + { + nElectronsVsMean->Write(); + } + if (diffVsDrift) + { + diffVsDrift->Write(); + } + if (diffPerSqrtLVsDrift) + { + diffPerSqrtLVsDrift->Write(); + } + if (diffDXVsDrift) + { + diffDXVsDrift->Write(); + } + if (diffDYVsDrift) + { + diffDYVsDrift->Write(); + } + if (diffDXPerSqrtLVsDrift) + { + diffDXPerSqrtLVsDrift->Write(); + } + if (diffDYPerSqrtLVsDrift) + { + diffDYPerSqrtLVsDrift->Write(); + } + if (driftXY) + { + driftXY->Write(); + } + if (electronDensityProfile) + { + electronDensityProfile->Write(); + } EDrift_outf->Close(); } + if (m_avgOutf && m_avgXResidualTree) + { + m_avgOutf->cd(); + m_avgXResidualTree->Write(); + if (Verbosity() > 0) + { + std::cout << Name() << ": wrote avgXResidual tree to " << m_avg_output_file << std::endl; + } + m_avgOutf->Close(); + } return Fun4AllReturnCodes::EVENT_OK; } @@ -936,29 +1731,75 @@ void PHG4TpcElectronDrift::set_seed(const unsigned int seed) void PHG4TpcElectronDrift::SetDefaultParameters() { - // longitudinal diffusion for 50:50 Ne:CF4 is 0.012, transverse is 0.004, drift velocity is 0.008 - // longitudinal diffusion for 60:40 Ar:CF4 is 0.012, transverse is 0.004, drift velocity is 0.008 (chosen to be the same at 50/50 Ne:CF4) - // longitudinal diffusion for 65:25:10 Ar:CF4:N2 is 0.013613, transverse is 0.005487, drift velocity is 0.006965 - // longitudinal diffusion for 75:20:05 Ar:CF4:i-C4H10 is 0.014596, transverse is 0.005313, drift velocity is 0.007550 + //longitudinal diffusion for 50:50 Ne:CF4 is 0.012, transverse is 0.004, drift velocity is 0.008 + //longitudinal diffusion for 60:40 Ar:CF4 is 0.012, transverse is 0.004, drift velocity is 0.008 (chosen to be the same at 50/50 Ne:CF4) + //longitudinal diffusion for 65:25:10 Ar:CF4:N2 is 0.013613, transverse is 0.005487, drift velocity is 0.006965 + //longitudinal diffusion for 75:20:05 Ar:CF4:i-C4H10 is 0.014596, transverse is 0.005313, drift velocity is 0.007550 set_default_double_param("diffusion_long", 0.014596); // cm/SQRT(cm) set_default_double_param("diffusion_trans", 0.005313); // cm/SQRT(cm) - set_default_double_param("Ne_frac", 0.00); - set_default_double_param("Ar_frac", 0.75); + set_default_double_param("drift_velocity", 0.00755); // cm/ns + set_default_double_param("Ne_frac", 0.00); + set_default_double_param("Ar_frac", 0.75); set_default_double_param("CF4_frac", 0.20); set_default_double_param("N2_frac", 0.00); set_default_double_param("isobutane_frac", 0.05); set_default_double_param("min_active_radius", 30.); // cm set_default_double_param("max_active_radius", 78.); // cm + set_default_double_param("max_time", 13200.); // ns + set_default_double_param("extended_readout_time", 7000.); // ns // These are purely fudge factors, used to increase the resolution to 150 microns and 500 microns, respectively // override them from the macro to get a different resolution set_default_double_param("added_smear_trans", 0.0); // cm (used to be 0.085 before sims got better) set_default_double_param("added_smear_long", 0.0); // cm (used to be 0.105 before sims got better) + set_default_double_param("force_min_trans_drift_length", 0.0); // cm + set_default_int_param("density_layer", -1); + set_default_double_param("density_bin_width_cm", 0.05); + set_default_double_param("density_window_cm", -1.0); + set_default_double_param("density_radial_margin_cm", 0.2); + set_default_int_param("uniform_density_test", 0); + set_default_int_param("average_x_layer", -1); return; } +bool PHG4TpcElectronDrift::cylinder_intersection(const double x0, const double y0, const double z0, + const double vx, const double vy, const double vz, + const double R, double &t_out, + double &xi, double &yi, double &zi) const +{ + const double a = square(vx) + square(vy); + const double b = 2.0 * (vx * x0 + vy * y0); + const double c = square(x0) + square(y0) - square(R); + const double disc = b * b - 4.0 * a * c; + if (a == 0.0 || disc < 0.0) + { + return false; + } + const double s = std::sqrt(disc); + const double t1 = (-b + s) / (2.0 * a); + const double t2 = (-b - s) / (2.0 * a); + double t = std::numeric_limits::max(); + if (t1 > 0.0) + { + t = std::min(t, t1); + } + if (t2 > 0.0) + { + t = std::min(t, t2); + } + if (!std::isfinite(t) || t <= 0.0) + { + return false; + } + t_out = t; + xi = x0 + t * vx; + yi = y0 + t * vy; + zi = z0 + t * vz; + return true; +} + void PHG4TpcElectronDrift::setTpcDistortion(PHG4TpcDistortion *distortionMap) { m_distortionMap.reset(distortionMap); @@ -983,7 +1824,7 @@ void PHG4TpcElectronDrift::registerPadPlane(PHG4TpcPadPlane *inpadplane) void PHG4TpcElectronDrift::set_flag_threshold_distortion(bool setflag, float setthreshold) { - std::cout << std::format("The logical status of threshold is now {}! and the value is set to {}", setflag, setthreshold) + std::cout << boost::str(boost::format("The logical status of threshold is now %d! and the value is set to %f") % setflag % setthreshold) << std::endl << std::endl << std::endl; diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h index de2810e918..aa21170be0 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h @@ -16,11 +16,14 @@ #include #include +#include #include #include #include #include #include +#include +#include class PHG4TpcPadPlane; class PHG4TpcDistortion; @@ -28,6 +31,7 @@ class PHCompositeNode; class TH1; class TH2; class TNtuple; +class TTree; class TFile; class TrkrHitSetContainer; class TrkrHitTruthAssoc; @@ -37,6 +41,7 @@ class TrkrTruthTrack; class DistortedTrackContainer; class TpcClusterBuilder; class PHG4TpcGeomContainer; +class SvtxTrackMap; class ClusHitsVerbose; class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface @@ -83,6 +88,15 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface void set_zero_bfield_flag(bool flag) { zero_bfield = flag; }; void set_zero_bfield_diffusion_factor(double f) { zero_bfield_diffusion_factor = f; }; void use_PDG_gas_params() { m_use_PDG_gas_params = true; } + void set_force_min_trans_drift_length(double length) + { + force_min_trans_drift_length = (length > 0.) ? length : 0.; + set_double_param("force_min_trans_drift_length", force_min_trans_drift_length); + } + void set_enable_laser_clustering(bool b) { m_enable_laser_clustering = b; } + void set_qa_output_file(const std::string &f) { m_qa_output_file = f; } + void set_avg_output_file(const std::string &f) { m_avg_output_file = f; } + void set_do_ElectronDriftQAHistos(bool b) { do_ElectronDriftQAHistos = b; } ClusHitsVerbosev1 *mClusHitsVerbose{nullptr}; private: @@ -93,17 +107,50 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface TrkrClusterContainer *truthclustercontainer{nullptr}; // the TrkrClusterContainer for truth clusters ActsGeometry *m_tGeometry{nullptr}; PHG4TpcGeomContainer *seggeo{nullptr}; + SvtxTrackMap *m_track_map{nullptr}; TNtuple *nt{nullptr}; TNtuple *nthit{nullptr}; TNtuple *ntfinalhit{nullptr}; TNtuple *ntpad{nullptr}; + TTree *m_avgXResidualTree{nullptr}; + float m_avgTree_trackid{0}; + float m_avgTree_residual{0}; + float m_avgTree_avgx{0}; + float m_avgTree_xint{0}; + float m_avgTree_residual_primary{0}; + float m_avgTree_avgx_primary{0}; + float m_avgTree_count_primary{0}; + float m_avgTree_residual_end{0}; + float m_avgTree_avgx_end{0}; + float m_avgTree_count_end{0}; + float m_avgTree_residual_end_primary{0}; + float m_avgTree_avgx_end_primary{0}; + float m_avgTree_count_end_primary{0}; + float m_avgTree_residual_x{0}; + float m_avgTree_avgx_cart{0}; + float m_avgTree_residual_x_primary{0}; + float m_avgTree_avgx_cart_primary{0}; + float m_avgTree_residual_x_end{0}; + float m_avgTree_avgx_cart_end{0}; + float m_avgTree_residual_x_end_primary{0}; + float m_avgTree_avgx_cart_end_primary{0}; ///@name evaluation histograms //@{ TH1 *dlong{nullptr}; TH1 *dtrans{nullptr}; TH1 *ratioElectronsRR{nullptr}; + TH1 *diffDistance{nullptr}; + TH1 *diffDX{nullptr}; + TH1 *diffDY{nullptr}; + TH1 *nElectrons{nullptr}; + TH1 *poissonMean{nullptr}; + TH1 *diffPerSqrtL{nullptr}; + TH1 *diffDXPerSqrtL{nullptr}; + TH1 *diffDYPerSqrtL{nullptr}; + TH1 *nElectronsPerCm{nullptr}; + TH1 *electronDensityProfile{nullptr}; TH2 *hitmapstart{nullptr}; TH2 *hitmapend{nullptr}; TH2 *hitmapstart_z{nullptr}; @@ -118,6 +165,14 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface TH2 *deltarnodiff{nullptr}; TH2 *deltarnodist{nullptr}; TH2 *deltaz{nullptr}; + TH2 *nElectronsVsMean{nullptr}; + TH2 *diffVsDrift{nullptr}; + TH2 *diffPerSqrtLVsDrift{nullptr}; + TH2 *diffDXVsDrift{nullptr}; + TH2 *diffDYVsDrift{nullptr}; + TH2 *diffDXPerSqrtLVsDrift{nullptr}; + TH2 *diffDYPerSqrtLVsDrift{nullptr}; + TNtuple *driftXY{nullptr}; //@} int event_num{0}; @@ -125,20 +180,69 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface float max_g4hitstep{7.}; float thresholdforreachesreadout{0.5}; - double diffusion_trans = std::numeric_limits::quiet_NaN(); - double added_smear_sigma_trans = std::numeric_limits::quiet_NaN(); - double diffusion_long = std::numeric_limits::quiet_NaN(); - double added_smear_sigma_long = std::numeric_limits::quiet_NaN(); - double tpc_length = std::numeric_limits::quiet_NaN(); - double electrons_per_gev = std::numeric_limits::quiet_NaN(); - double min_active_radius = std::numeric_limits::quiet_NaN(); - double max_active_radius = std::numeric_limits::quiet_NaN(); - double min_time = std::numeric_limits::quiet_NaN(); - double max_time = std::numeric_limits::quiet_NaN(); + double diffusion_trans = std::numeric_limits::signaling_NaN(); + double added_smear_sigma_trans = std::numeric_limits::signaling_NaN(); + double diffusion_long = std::numeric_limits::signaling_NaN(); + double added_smear_sigma_long = std::numeric_limits::signaling_NaN(); + double drift_velocity = std::numeric_limits::signaling_NaN(); + double tpc_length = std::numeric_limits::signaling_NaN(); + double electrons_per_gev = std::numeric_limits::signaling_NaN(); + double min_active_radius = std::numeric_limits::signaling_NaN(); + double max_active_radius = std::numeric_limits::signaling_NaN(); + double min_time = std::numeric_limits::signaling_NaN(); + double max_time = std::numeric_limits::signaling_NaN(); double zero_bfield_diffusion_factor{3.5}; // at drift voltage of 400 V + double force_min_trans_drift_length{0.0}; + bool m_density_enabled{false}; + int m_density_layer{-1}; + double m_density_bin_width_cm{0.05}; + double m_density_window_cm{0.5}; + double m_density_layer_radius{0.0}; + double m_density_layer_rlow{0.0}; + double m_density_layer_rhigh{0.0}; + double m_density_radial_margin_cm{0.0}; + double m_density_capture_rlow{0.0}; + double m_density_capture_rhigh{0.0}; + double m_density_hist_half_range{0.0}; + bool m_avg_x_enabled{false}; + int m_avg_layer{-1}; + double m_avg_layer_radius{0.0}; + double m_avg_layer_rlow{0.0}; + double m_avg_layer_rhigh{0.0}; + std::unordered_map m_track_path_offset; + std::unordered_map m_track_anchor_set; + std::unordered_map m_track_anchor_length; + struct TrackLayerData + { + double sum_rphi{0.0}; + double sum_rphi_primary{0.0}; + double sum_rphi_end{0.0}; + double sum_rphi_end_primary{0.0}; + double sum_x{0.0}; + double sum_x_primary{0.0}; + double sum_x_end{0.0}; + double sum_x_end_primary{0.0}; + std::size_t count{0}; + std::size_t count_primary{0}; + std::size_t count_end{0}; + std::size_t count_end_primary{0}; + double base_x{0.0}; + double base_y{0.0}; + double base_z{0.0}; + double dir_x{0.0}; + double dir_y{0.0}; + double dir_z{0.0}; + bool have_line{false}; + }; + std::unordered_map m_track_layer_data; + bool m_uniform_density_test{false}; + std::vector cluster_size_cdf; + bool m_enable_laser_clustering{false}; + std::string m_qa_output_file{"ElectronDriftQA.root"}; + std::string m_avg_output_file{"avgXResidual.root"}; bool record_ClusHitsVerbose{false}; - bool do_ElectronDriftQAHistos{false}; + bool do_ElectronDriftQAHistos{true}; bool do_getReachReadout{false}; bool zero_bfield{false}; bool m_use_PDG_gas_params{false}; @@ -149,6 +253,7 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface std::unique_ptr m_distortionMap; std::unique_ptr m_outf; std::unique_ptr EDrift_outf; + std::unique_ptr m_avgOutf; std::string detector; std::string hitnodename; @@ -162,6 +267,11 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface void operator()(gsl_rng *rng) const { gsl_rng_free(rng); } }; std::unique_ptr RandomGenerator; + + bool cylinder_intersection(double x0, double y0, double z0, + double vx, double vy, double vz, + double R, double &t_out, + double &xi, double &yi, double &zi) const; }; #endif // G4TPC_PHG4TPCELECTRONDRIFT_H diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index 431f9e684c..f0060ee64d 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -10,10 +10,8 @@ #include #include - #include #include - // Move to new storage containers #include #include // for hitkey, hitse... @@ -29,12 +27,9 @@ #include // for PHWHERE -#include #include #include -#include - - +#include #include #include #include @@ -42,17 +37,22 @@ #include #include #include - +#include #include #include // for gsl_rng_alloc +#include + +#include +#include #include #include // for getenv #include #include #include // for _Rb_tree_cons... #include // for pair +#include // for std::ifstream #include #include @@ -63,7 +63,7 @@ namespace { //! convenient square function template - constexpr T square(const T &x) + inline constexpr T square(const T &x) { return x * x; } @@ -81,30 +81,67 @@ namespace return std::exp(-square(x / sigma) / 2) / (sigma * std::sqrt(2 * M_PI)); } - constexpr unsigned int print_layer = 18; + static constexpr int print_layer = 18; } // namespace PHG4TpcPadPlaneReadout::PHG4TpcPadPlaneReadout(const std::string &name) : PHG4TpcPadPlane(name) - , RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937)) { InitializeParameters(); // if(m_flagToUseGain==1) ReadGain(); - + RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937); gsl_rng_set(RandomGenerator, PHRandomSeed()); // fixed seed is handled in this funtcion + return; } PHG4TpcPadPlaneReadout::~PHG4TpcPadPlaneReadout() { gsl_rng_free(RandomGenerator); - for (auto *his : h_gain) + for (auto his : h_gain) { delete his; } + + for(int i=0; i<2; ++i) + for(int j=0; j<3; ++j) + for(int k=0; k<12; ++k) + delete flangau[i][j][k]; +} + +// pick all layers whose radial annulus intersects [rad - nsig*sigma, rad + nsig*sigma] +std::vector +PHG4TpcPadPlaneReadout::layersInRadialWindow(double rad, double sigma, double nsig) const +{ + std::vector out; + const double rmin = rad - nsig*sigma; + const double rmax = rad + nsig*sigma; + + PHG4TpcGeomContainer::ConstRange layerrange = GeomContainer->get_begin_end(); + for (auto it = layerrange.first; it != layerrange.second; ++it) + { + const auto* g = it->second; + const double rLow = g->get_radius() - 0.5*g->get_thickness(); + const double rHigh = g->get_radius() + 0.5*g->get_thickness(); + + const double lo = std::max(rLow, rmin); + const double hi = std::min(rHigh, rmax); + if (hi > lo) out.push_back(g->get_layer()); + } + return out; +} + +// get geometry for a given layer (utility; linear scan is fine here) +PHG4TpcGeom* +PHG4TpcPadPlaneReadout::getGeomForLayer(unsigned int layer) const +{ + PHG4TpcGeomContainer::ConstRange layerrange = GeomContainer->get_begin_end(); + for (auto it = layerrange.first; it != layerrange.second; ++it) + if (static_cast(it->second->get_layer()) == layer) return it->second; + return nullptr; } //_________________________________________________________ @@ -119,71 +156,53 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) const std::string seggeonodename = "TPCGEOMCONTAINER"; GeomContainer = findNode::getClass(topNode, seggeonodename); assert(GeomContainer); - - PHG4TpcGeom *layergeom = GeomContainer->GetLayerCellGeom(20); // z geometry is the same for all layers - double tpc_adc_clock = layergeom->get_adc_clock(); - double extended_readout_time = layergeom->get_extended_readout_time(); - double maxdriftlength = layergeom->get_max_driftlength(); - double drift_velocity_sim = layergeom->get_drift_velocity_sim(); - const double TBinWidth = tpc_adc_clock; - const double MaxT = extended_readout_time + 2.0 * maxdriftlength / drift_velocity_sim; // allows for extended time readout - const double MinT = 0; - NTBins = (int) ((MaxT - MinT) / TBinWidth) + 1; - - if (m_use_module_gain_weights) - { - int side; - int region; - int sector; - double weight; - std::ifstream weights_file(m_tpc_module_gain_weights_file); - if (!weights_file.is_open()) + if(m_use_module_gain_weights) { - std::cout << ".In PHG4TpcPadPlaneReadout: Option to use module gain weights enabled, but weights file not found. Aborting." << std::endl; - return Fun4AllReturnCodes::ABORTEVENT; - } + int side, region, sector; + double weight; + std::ifstream weights_file(m_tpc_module_gain_weights_file); + if(!weights_file.is_open()) + { + std::cout << ".In PHG4TpcPadPlaneReadout: Option to use module gain weights enabled, but weights file not found. Aborting." << std::endl; + return Fun4AllReturnCodes::ABORTEVENT; + } - for (int iside = 0; iside < 2; ++iside) - { - for (int isec = 0; isec < 12; ++isec) - { - for (int ir = 0; ir < 3; ++ir) - { - weights_file >> side >> region >> sector >> weight; - m_module_gain_weight[side][region][sector] = weight; - std::cout << " iside " << iside << " side " << side << " ir " << ir - << " region " << region << " isec " << isec - << " sector " << sector << " weight " << weight << std::endl; - } - } - } - } + for(int iside =0; iside < 2; ++iside) + { + for(int isec = 0; isec < 12; ++isec) + { + for(int ir = 0; ir < 3; ++ir) + { + weights_file >> side >> region >> sector >> weight; + m_module_gain_weight[side][region][sector] = weight; + std::cout << " iside " << iside << " side " << side << " ir " << ir + << " region " << region << " isec " << isec + << " sector " << sector << " weight " << weight << std::endl; + } + } + } + } - if (m_useLangau) - { - int side; - int region; - int sector; - double par0; - double par1; - double par2; - double par3; - std::ifstream pars_file(m_tpc_langau_pars_file); - if (!pars_file.is_open()) + if(m_useLangau) { - std::cout << ".In PHG4TpcPadPlaneReadout: Option to use Langau parameters enabled, but parameter file not found. Aborting." << std::endl; - return Fun4AllReturnCodes::ABORTEVENT; - } + int side, region, sector; + double par0; double par1; double par2; double par3; + std::ifstream pars_file(m_tpc_langau_pars_file); + if(!pars_file.is_open()) + { + std::cout << ".In PHG4TpcPadPlaneReadout: Option to use Langau parameters enabled, but parameter file not found. Aborting." << std::endl; + return Fun4AllReturnCodes::ABORTEVENT; + } - for (int iside = 0; iside < 2; ++iside) - { - for (int isec = 0; isec < 12; ++isec) - { - for (int ir = 0; ir < 3; ++ir) - { - pars_file >> side >> region >> sector >> par0 >> par1 >> par2 >> par3; - flangau[side][region][sector] = new TF1(std::format("flangau_{}_{}_{}", side, region, sector).c_str(), [](double *x, double *par) - { + for(int iside =0; iside < 2; ++iside) + { + for(int isec = 0; isec < 12; ++isec) + { + for(int ir = 0; ir < 3; ++ir) + { + pars_file >> side >> region >> sector >> par0 >> par1 >> par2 >> par3; + flangau[side][region][sector] = new TF1((boost::format("flangau_%d_%d_%d") % side % region % sector).str().c_str(), [](double *x, double *par) + { Double_t invsq2pi = 0.3989422804014; Double_t mpshift = -0.22278298; Double_t np = 100.0; @@ -192,8 +211,7 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) Double_t mpc; Double_t fland; Double_t sum = 0.0; - Double_t xlow; - Double_t xupp; + Double_t xlow,xupp; Double_t step; Double_t i; mpc = par[1] - mpshift * par[0]; @@ -211,17 +229,19 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) sum += fland * TMath::Gaus(x[0],xx,par[3]); } - return (par[2] * step * sum * invsq2pi / par[3]); }, 0, 5000, 4); + return (par[2] * step * sum * invsq2pi / par[3]); + }, 0, 5000, 4); - flangau[side][region][sector]->SetParameters(par0, par1, par2, par3); - // std::cout << " iside " << iside << " side " << side << " ir " << ir - // << " region " << region << " isec " << isec - // << " sector " << sector << " weight " << weight << std::endl; - } - } + + flangau[side][region][sector]->SetParameters(par0,par1,par2,par3); + //std::cout << " iside " << iside << " side " << side << " ir " << ir + // << " region " << region << " isec " << isec + // << " sector " << sector << " weight " << weight << std::endl; + } + } + } } - } - if (m_maskDeadChannels) + if (m_maskDeadChannels) { makeChannelMask(m_deadChannelMap, m_deadChannelMapName, "TotalDeadChannels"); } @@ -229,22 +249,102 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) { makeChannelMask(m_hotChannelMap, m_hotChannelMapName, "TotalHotChannels"); } - - std::cout<<"!!!!!Before Load Maps"< missing_layers; + + // We consider readout layers to start at 7 (Pads is sized as 7 + 16*3) + constexpr unsigned int kFirstPolygonLayer = 7; - for (int j=0; j<7+16*3;j++){ - for(size_t i = 0; i < Pads[j].size(); i++) + // Check if any layer has polygons at all + for (const auto &layerPads : Pads) { - std::cout<<"Module "<<(j-7)/16<<" layer = "<= 3) { poly_layers++; break; } + } + } + m_serf_polygons_present = (poly_layers > 0); + + // If SERF is desired, verify full coverage of all readout layers present in geometry + if (m_use_serf_padsharing) + { + if (!GeomContainer) + { + std::cerr << "PHG4TpcPadPlaneReadout: Geometry container missing in InitRun." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + PHG4TpcGeomContainer::ConstRange layerrange = GeomContainer->get_begin_end(); + for (auto it = layerrange.first; it != layerrange.second; ++it) + { + unsigned int layer = static_cast(it->second->get_layer()); + if (layer < kFirstPolygonLayer) continue; // allow early layers to miss polygons + bool has_poly = false; + if (layer < Pads.size()) + { + const auto &layerPads = Pads[layer]; + for (const auto &pad : layerPads) + { + if (pad.vertices.size() >= 3) { has_poly = true; break; } + } + } + if (!has_poly) missing_layers.push_back(layer); + } + + if (!m_serf_polygons_present) + { + std::cout << "PHG4TpcPadPlaneReadout: SERF requested but pad polygons not available; using analytic sharing." + << " Call RequireSerfPadSharing(true) to abort instead." << std::endl; + if (m_require_serf) + { + std::cerr << "PHG4TpcPadPlaneReadout: SERF pad polygons required but missing. Aborting run." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + } + + if (!missing_layers.empty()) + { + std::cout << "PHG4TpcPadPlaneReadout: SERF polygons missing for " + << missing_layers.size() << " readout layer(s) starting at layer 7." << std::endl; + if (m_require_serf_full) + { + std::cerr << "PHG4TpcPadPlaneReadout: RequireSerfFullCoverage enabled; aborting due to missing polygons in readout layers." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + } } } + // 3) (optional) print a summary + + + + + /* for (int j=0; j<7+16*3;j++){ + for(size_t i = 0; i < Pads[j].size(); i++) + { + std::cout<<"Module "<<(j-7)/16<<" layer = "< + PHG4TpcPadPlaneReadout::brdMaps_ = { + "/sphenix/user/mitrankova/Simulation/PadPlane/AutoPad-R1-RevA.brd", + "/sphenix/user/mitrankova/Simulation/PadPlane/AutoPad-R2-RevA-Pads.brd", + "/sphenix/user/mitrankova/Simulation/PadPlane/AutoPad-R3-RevA.brd" +}; + //_________________________________________________________ void PHG4TpcPadPlaneReadout::loadPadPlanes() { @@ -354,7 +454,6 @@ void PHG4TpcPadPlaneReadout::loadPadPlanes() { } } - //_________________________________________________________ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification() { @@ -368,11 +467,11 @@ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification() // and yes, the parameter you're looking for is of course the slope, which is the inverse gain. double nelec = gsl_ran_exponential(RandomGenerator, averageGEMGain); if (m_usePolya) - { + { double y; double xmax = 5000; double ymax = 0.376; - while (true) + while (true) { nelec = gsl_ran_flat(RandomGenerator, 0, xmax); y = gsl_rng_uniform(RandomGenerator) * ymax; @@ -405,11 +504,11 @@ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification(double weight) double y; double xmax = 5000; double ymax = 0.376; - while (true) + while (true) { nelec = gsl_ran_flat(RandomGenerator, 0, xmax); y = gsl_rng_uniform(RandomGenerator) * ymax; - if (y <= pow((1 + polyaTheta) * (nelec / q_bar), polyaTheta) * exp(-(1 + polyaTheta) * (nelec / q_bar))) + if (y <= pow((1 + polyaTheta) * (nelec / q_bar), polyaTheta) * exp(-(1 + polyaTheta) * (nelec / q_bar))) { break; } @@ -423,41 +522,130 @@ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification(double weight) //_________________________________________________________ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification(TF1 *f) { - double nelec = f->GetRandom(0, 5000); + double nelec = f->GetRandom(0,5000); // Put gain reading here return nelec; } -//_________________________________________________________ -inline void rotatePointToSector(double x, double y, - double& xNew, double& yNew, - int& sector) + + +void PHG4TpcPadPlaneReadout::rotatePointToSector( + double x, double y, + unsigned int side, + int& sectorFound, + double& xNew, double& yNew +) { - // 1) compute original phi in [−π, +π] - double phi = std::atan2(y, x); - - // 2) find the 30°‐wide wedge it lives in - const double PI = std::acos(-1.0); - double wedgeWidth = 2.0 * PI / TpcDefs::NSectors; // = π/6 - // shift by half‐wedge so floor() bins correctly - sector = static_cast( - std::floor((phi + wedgeWidth * 0.5) / wedgeWidth) - ) % TpcDefs::NSectors; - if (sector < 0) sector += TpcDefs::NSectors; // ensure non‐negative - - // 3) how much to rotate so that this sector’s center → +90° (π/2) - double targetCenter = PI / 2.0; // 12 o'clock - double originalCenter = sector * wedgeWidth; // e.g. 3 → π/2 - double dphi = targetCenter - originalCenter; - - // 4) apply rotation in polar coords - double R = std::hypot(x, y); - double phiRot = phi + dphi; - xNew = R * std::cos(phiRot); - yNew = R * std::sin(phiRot); + // ---- helpers ------------------------------------------------------------- + const double PI = std::acos(-1.0); + const double TWOPI = 2.0*PI; + + auto wrap = [&](double a) { + while (a <= -PI) a += TWOPI; + while (a > PI) a -= TWOPI; + return a; + }; + + // check a ∈ [lo,hi) with wrap at ±π + auto inInterval = [&](double a, double lo, double hi) { + a = wrap(a); lo = wrap(lo); hi = wrap(hi); + if (lo <= hi) return (a >= lo && a < hi); + // interval crosses the branch cut + return (a >= lo || a < hi); + }; + + // midpoint on the circle between lo..hi (shorter arc) + auto mid = [&](double lo, double hi) { + double d = wrap(hi - lo); + return wrap(lo + 0.5*d); + }; + + // ---- 1) angle & sector --------------------------------------------------- + const double R = std::hypot(x, y); + const double phi = wrap(std::atan2(y, x)); + + // sector_min_Phi / sector_max_Phi are class members filled from LayerGeom + sectorFound = -1; + for (int s = 0; s < 12; ++s) + { + if (inInterval(phi, sector_min_Phi[side][s], sector_max_Phi[side][s])) + { sectorFound = s; break; } + } + + // If exactly on a boundary, pick nearest sector center + if (sectorFound < 0) + { + double best = 1e9; int bestS = 0; + for (int s = 0; s < 12; ++s) + { + double c = mid(sector_min_Phi[side][s], sector_max_Phi[side][s]); + double d = std::fabs(wrap(phi - c)); + if (d < best) { best = d; bestS = s; } + } + sectorFound = bestS; + } + + // ---- 2) rotate to the reference sector (sector 2) ------------------------ + // The SERF reference polygons are defined in the frame of sector 2. + // Rotate points from the found sector into the sector-2 frame by the + // boundary difference: + // side 0 (South): dphi = sec_min[found] - sec_min[2] + // side 1 (North): dphi = sec_max[found] - sec_max[2] + // We then rotate the point by -dphi so that boundaries align. + const int refSector = 2; + const double currentBoundary = (side == 0) + ? sector_min_Phi[side][sectorFound] + : sector_max_Phi[side][sectorFound]; + const double refBoundary = (side == 0) + ? sector_min_Phi[side][refSector] + : sector_max_Phi[side][refSector]; + const double dphi = wrap(currentBoundary - refBoundary); + const double phiRot = wrap(phi - dphi); + + xNew = R * std::cos(phiRot); + yNew = R * std::sin(phiRot); + + // ---- 3) mirror for South side to match reference polygon handedness ------ + // Looking from the South side outward (toward +z), the x-axis is to the left + // while the reference sector has x to the right, so flip x for side 0. + constexpr unsigned SOUTH_SIDE = 0; + if (side == SOUTH_SIDE) + { + xNew = -xNew; + } } + + + + +/* //_________________________________________________________ +std::vector computeShaperKernel() { + int NT = static_cast(DetectorParams::window_ns / DetectorParams::adc_dt); + std::vector h(NT); + for(int i = 0; i < NT; ++i) { + double t = (i + 0.5) * DetectorParams::adc_dt; + h[i] = (t / std::pow(DetectorParams::tau_shaper,2)) * std::exp(-t / DetectorParams::tau_shaper); + } + // normalize + double sum = 0; + for(double v : h) sum += v * DetectorParams::adc_dt; + for(double &v : h) v /= sum; + return h; +} + +//_________________________________________________________ + +double gaussianIntegral1D(double a, double b, double mu, double sigma) { + if(sigma <= 0) return 0.0; + return 0.5 * (TMath::Erf((b-mu)/(std::sqrt(2)*sigma)) - TMath::Erf((a-mu)/(std::sqrt(2)*sigma))); +} + +//_________________________________________________________ + +*/ + bool PHG4TpcPadPlaneReadout::pointInPolygon( double x, double y, const std::vector& poly @@ -473,8 +661,6 @@ bool PHG4TpcPadPlaneReadout::pointInPolygon( // does edge (j→i) straddle the horizontal ray at y? bool cond = ((yi > y) != (yj > y)); if (cond) { - // std::cout<<"Checking point ("<= 15 && (phi_gain * 180.0 / M_PI) < 345) + + if( (phi_gain*180.0/M_PI) >=15 && (phi_gain*180.0 / M_PI) < 345) { - sector = 1 + (int) ((phi_gain * 180.0 / M_PI - 15) / phistep); - } + sector = 1 + (int) ( (phi_gain*180.0/M_PI - 15) / phistep); + } else { sector = 0; } - + int this_region = -1; for (int iregion = 0; iregion < 3; ++iregion) { if (rad_gem < MaxRadius[iregion] && rad_gem > MinRadius[iregion]) { - this_region = iregion; + this_region = iregion; } } - if (this_region > -1) + if(this_region > -1) { nelec = getSingleEGEMAmplification(flangau[side][this_region][sector]); } - else + else { nelec = getSingleEGEMAmplification(); } } - + // std::cout<<"PHG4TpcPadPlaneReadout::MapToPadPlane gain_weight = "< pad_layer; - std::vector pad_phibin; - std::vector pad_phibin_share; - // populate_zigzag_phibins(side, layernum, phi, sigmaT, pad_phibin, pad_phibin_share); - SERF_zigzag_phibins(side, layernum, phi, rad_gem, sigmaT, pad_layer, pad_phibin, pad_phibin_share); - /* if (pad_phibin.size() == 0) { */ - /* pass_data.neff_electrons = 0; */ - /* } else { */ - /* pass_data.fillPhiBins(pad_phibin); */ - /* } */ - // Normalize the shares so they add up to 1 - double norm1 = 0.0; - for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad) + // NOTE: we no longer compute single-layer pad sharing here. Instead, we + // compute per-layer pad shares inside the candidate-layers loop below + // (using SERF_zigzag_phibins) to avoid double counting and to include + // all pads within 5σ radially across layers. +/* +norm1 = 0.0; + for (unsigned int ipad = 0; ipad < pad_phibin_ref.size(); ++ipad) { - double pad_share = pad_phibin_share[ipad]; + double pad_share = pad_phibin_share_ref[ipad]; norm1 += pad_share; } - for (unsigned int iphi = 0; iphi < pad_phibin.size(); ++iphi) + for (unsigned int iphi = 0; iphi < pad_phibin_ref.size(); ++iphi) { - pad_phibin_share[iphi] /= norm1; - } + pad_phibin_share_ref[iphi] /= norm1; + }*/ + /*std::cout<<"--------------------------"< 100 && layernum == print_layer) - { - std::cout << " populate t bins for layernum " << layernum - << " with t_gem " << t_gem << " SAMPA peaking time " << Ts << std::endl; - } + { + std::cout << " populate t bins for layernum " << layernum + << " with t_gem " << t_gem << " sigmaL[0] " << sigmaL[0] << " sigmaL[1] " << sigmaL[1] << std::endl; + } std::vector adc_tbin; std::vector adc_tbin_share; sampaTimeDistribution(t_gem, adc_tbin, adc_tbin_share); - /* if (adc_tbin.size() == 0) { */ /* pass_data.neff_electrons = 0; */ /* } else { */ @@ -742,141 +954,250 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( { adc_tbin_share[it] /= tnorm; } - - /* - if(layernum == print_layer) - { - std::cout << "t_gem " << t_gem << std::endl; - for (unsigned int it = 0; it < adc_tbin.size(); ++it) - { - std::cout << " tbin " << adc_tbin[it] << " share " << adc_tbin_share[it] << std::endl; - } - } - */ - - // Fill HitSetContainer - //=============== - // These are used to do a quick clustering for checking + // accumulate quick centroid info across all candidate layers double phi_integral = 0.0; double t_integral = 0.0; double weight = 0.0; - for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad) + // Determine if SERF pad sharing is possible (pad polygons loaded) + bool have_polygons = true; + for (unsigned int layer_cand : cand_layers) { - int pad_num = pad_phibin[ipad]; - double pad_share = pad_phibin_share[ipad]; - const unsigned int layer_i = pad_layer[ipad]; - - for (unsigned int it = 0; it < adc_tbin.size(); ++it) + if (layer_cand >= Pads.size() || Pads[layer_cand].empty()) { - int tbin_num = adc_tbin[it]; - double adc_bin_share = adc_tbin_share[it]; + have_polygons = false; + break; + } + } - // Divide electrons from avalanche between bins - float neffelectrons = nelec * (pad_share) * (adc_bin_share); - if (neffelectrons < neffelectrons_threshold) - { - continue; // skip signals that will be below the noise suppression threshold - } + if (m_use_serf_padsharing && have_polygons) + { + // Pass 1: compute total Gaussian mass across all candidate layers + double total_mass = 0.0; + for (unsigned int layer_cand : cand_layers) + { + PHG4TpcGeom* thisGeom = getGeomForLayer(layer_cand); + if (!thisGeom) continue; + LayerGeom = thisGeom; + sector_min_Phi = LayerGeom->get_sector_min_phi(); + sector_max_Phi = LayerGeom->get_sector_max_phi(); + phi_bin_width = LayerGeom->get_phistep(); + const double phi_for_layer = check_phi(side, phi, rad_gem); + + std::vector pad_phibin_tmp; + std::vector pad_mass_tmp; // unnormalized integrated mass per pad + SERF_zigzag_phibins(side, layer_cand, phi_for_layer, rad_gem, sigmaT, pad_phibin_tmp, pad_mass_tmp); + for (double v : pad_mass_tmp) total_mass += v; + } - if (tbin_num >= tbins) - { - std::cout << " Error making key: adc_tbin " << tbin_num << " ntbins " << tbins << std::endl; - } - if (pad_num >= phibins) - { - std::cout << " Error making key: pad_phibin " << pad_num << " nphibins " << phibins << std::endl; - } + if (total_mass <= 1e-16) total_mass = 1.0; // avoid division by zero - // collect information to do simple clustering. Checks operation of PHG4CylinderCellTpcReco, and - // is also useful for comparison with PHG4TpcClusterizer result when running single track events. - // The only information written to the cell other than neffelectrons is tbin and pad number, so get those from geometry - double tcenter = LayerGeom->get_zcenter(tbin_num); - double phicenter = LayerGeom->get_phicenter(pad_num, side); - phi_integral += phicenter * neffelectrons; - t_integral += tcenter * neffelectrons; - weight += neffelectrons; - if (Verbosity() > 1 && layernum == print_layer) + // Pass 2: fill hits with correct per-layer mass fraction + for (unsigned int layer_cand : cand_layers) + { + PHG4TpcGeom* thisGeom = getGeomForLayer(layer_cand); + if (!thisGeom) continue; + LayerGeom = thisGeom; + sector_min_Phi = LayerGeom->get_sector_min_phi(); + sector_max_Phi = LayerGeom->get_sector_max_phi(); + phi_bin_width = LayerGeom->get_phistep(); + const double phi_for_layer = check_phi(side, phi, rad_gem); + + std::vector pad_phibin; + std::vector pad_mass; // unnormalized integrated mass per pad + SERF_zigzag_phibins(side, layer_cand, phi_for_layer, rad_gem, sigmaT, pad_phibin, pad_mass); + + const auto phibins = LayerGeom->get_phibins(); + const auto tbins = LayerGeom->get_zbins(); + const unsigned int pads_per_sector = phibins / 12; + + for (size_t i = 0; i < pad_phibin.size(); ++i) { - std::cout << " tbin_num " << tbin_num << " tcenter " << tcenter << " pad_num " << pad_num << " phicenter " << phicenter - << " neffelectrons " << neffelectrons << " neffelectrons_threshold " << neffelectrons_threshold << std::endl; - } + const int pad_num = pad_phibin[i]; + const double pad_fraction = pad_mass[i] / total_mass; // fraction of total mass + if (pad_fraction <= 0) continue; - // new containers - //============ - // We add the Tpc TrkrHitsets directly to the node using hitsetcontainer - // We need to create the TrkrHitSet if not already made - each TrkrHitSet should correspond to a Tpc readout module - // The hitset key includes the layer, sector, side + const unsigned int sector = (pad_num >= 0) ? (static_cast(pad_num) / pads_per_sector) : 0; + TrkrDefs::hitsetkey hitsetkey = TpcDefs::genHitSetKey(layer_cand, sector, side); + auto hitsetit = hitsetcontainer->findOrAddHitSet(hitsetkey); + auto single_hitsetit = single_hitsetcontainer->findOrAddHitSet(hitsetkey); - // The side is an input parameter + for (unsigned int itb = 0; itb < adc_tbin.size(); ++itb) + { + const int tbin_num = adc_tbin[itb]; + if (tbin_num < 0 || tbin_num >= static_cast(tbins)) continue; + if (pad_num < 0 || pad_num >= static_cast(phibins)) continue; + + const double tshare = adc_tbin_share[itb]; + const float neffelectrons_bin = nelec * pad_fraction * tshare; + if (neffelectrons_bin < neffelectrons_threshold) continue; + + TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(static_cast(pad_num), static_cast(tbin_num)); + TrkrHit* hit = hitsetit->second->getHit(hitkey); + if (!hit) { hit = new TrkrHitv2(); hitsetit->second->addHitSpecificKey(hitkey, hit); } + hit->addEnergy(neffelectrons_bin); + + TrkrHit* single_hit = single_hitsetit->second->getHit(hitkey); + if (!single_hit) { single_hit = new TrkrHitv2(); single_hitsetit->second->addHitSpecificKey(hitkey, single_hit); } + single_hit->addEnergy(neffelectrons_bin); + + if (Verbosity() > 50) + { + const int layer_print = static_cast(layer_cand); + if (m_dbg_pad_hit_layer < 0 || layer_print == m_dbg_pad_hit_layer) + { + std::cout << "PadHit: layer=" << layer_cand + << " side=" << side + << " pad=" << pad_num + << " tbin=" << tbin_num + << " energy=" << neffelectrons_bin + << std::endl; + } + } - // get the Tpc readout sector - there are 12 sectors with how many pads each? - unsigned int pads_per_sector = phibins / 12; - unsigned int sector = pad_num / pads_per_sector; - TrkrDefs::hitsetkey hitsetkey = TpcDefs::genHitSetKey(layer_i, sector, side); - // Use existing hitset or add new one if needed - TrkrHitSetContainer::Iterator hitsetit = hitsetcontainer->findOrAddHitSet(hitsetkey); - TrkrHitSetContainer::Iterator single_hitsetit = single_hitsetcontainer->findOrAddHitSet(hitsetkey); - TrkrDefs::hitkey hitkey; + tpc_truth_clusterer.addhitset(hitsetkey, hitkey, neffelectrons_bin); - if (m_maskDeadChannels) - { - hitkey = TpcDefs::genHitKey((unsigned int) pad_num, 0); - if (m_deadChannelMap.contains(hitsetkey) && - std::find(m_deadChannelMap[hitsetkey].begin(), m_deadChannelMap[hitsetkey].end(), hitkey) != m_deadChannelMap[hitsetkey].end()) - { - continue; + const double tcenter = LayerGeom->get_zcenter(tbin_num); + const double phicenter = LayerGeom->get_phicenter(pad_num, side); + phi_integral += phicenter * neffelectrons_bin; + t_integral += tcenter * neffelectrons_bin; + weight += neffelectrons_bin; } } - if (m_maskHotChannels) - { - hitkey = TpcDefs::genHitKey((unsigned int) pad_num, 0); - if (m_hotChannelMap.contains(hitsetkey) && - std::find(m_hotChannelMap[hitsetkey].begin(), m_hotChannelMap[hitsetkey].end(), hitkey) != m_hotChannelMap[hitsetkey].end()) - { - continue; - } - } - // generate the key for this hit, requires tbin and phibin - hitkey = TpcDefs::genHitKey((unsigned int) pad_num, (unsigned int) tbin_num); + } + } + else + { + if (m_use_serf_padsharing && !have_polygons && !m_warned_serf_fallback) + { + std::cout << "PHG4TpcPadPlaneReadout: SERF polygons incomplete for some layers; per-hit fallback to analytic sharing is active." << std::endl; + m_warned_serf_fallback = true; + } + // Analytic triangular response across all layers within 5σ + // Pass 1: compute total radial weight across candidate layers + const double inv_sqrt2_sigma = 1.0 / (M_SQRT2 * sigmaT); + double total_radw = 0.0; + for (unsigned int layer_cand : cand_layers) + { + PHG4TpcGeom* thisGeom = getGeomForLayer(layer_cand); + if (!thisGeom) continue; + const double rcen = thisGeom->get_radius(); + const double thk = thisGeom->get_thickness(); + const double rlow = rcen - 0.5*thk; + const double rhigh= rcen + 0.5*thk; + const double dx1 = (rhigh - rad_gem) * inv_sqrt2_sigma; + const double dx0 = (rlow - rad_gem) * inv_sqrt2_sigma; + const double radw = 0.5 * (std::erf(dx1) - std::erf(dx0)); + if (radw > 0) total_radw += radw; + } + if (total_radw <= 1e-16) total_radw = 1.0; - // See if this hit already exists - TrkrHit *hit = nullptr; - hit = hitsetit->second->getHit(hitkey); - if (!hit) + // Pass 2: per-layer phi sharing and fill + for (unsigned int layer_cand : cand_layers) + { + PHG4TpcGeom* thisGeom = getGeomForLayer(layer_cand); + if (!thisGeom) continue; + LayerGeom = thisGeom; + sector_min_Phi = LayerGeom->get_sector_min_phi(); + sector_max_Phi = LayerGeom->get_sector_max_phi(); + phi_bin_width = LayerGeom->get_phistep(); + + const double rcen = LayerGeom->get_radius(); + const double thk = LayerGeom->get_thickness(); + const double rlow = rcen - 0.5*thk; + const double rhigh= rcen + 0.5*thk; + const double dx1 = (rhigh - rad_gem) * inv_sqrt2_sigma; + const double dx0 = (rlow - rad_gem) * inv_sqrt2_sigma; + const double radw = std::max(0.0, 0.5 * (std::erf(dx1) - std::erf(dx0))); + if (radw <= 0) continue; + + const double phi_for_layer = check_phi(side, phi, rad_gem); + std::vector pad_phibin; + std::vector pad_share_phi; + populate_zigzag_phibins(side, layer_cand, phi_for_layer, sigmaT, pad_phibin, pad_share_phi); + + double norm_phi = 0.0; for (double v : pad_share_phi) norm_phi += v; + if (norm_phi <= 1e-16) continue; + + const auto phibins = LayerGeom->get_phibins(); + const auto tbins = LayerGeom->get_zbins(); + const unsigned int pads_per_sector = phibins / 12; + + for (size_t i = 0; i < pad_phibin.size(); ++i) { - // create a new one - hit = new TrkrHitv2(); - hitsetit->second->addHitSpecificKey(hitkey, hit); - } - // Either way, add the energy to it -- adc values will be added at digitization - hit->addEnergy(neffelectrons); + const int pad_num = pad_phibin[i]; + const double pshare_phi = pad_share_phi[i] / norm_phi; + const double pad_fraction = (radw / total_radw) * pshare_phi; + if (pad_fraction <= 0) continue; + + const unsigned int sector = (pad_num >= 0) ? (static_cast(pad_num) / pads_per_sector) : 0; + TrkrDefs::hitsetkey hitsetkey = TpcDefs::genHitSetKey(layer_cand, sector, side); + if (m_maskDeadChannels) +{ + TrkrDefs::hitkey checkkey = TpcDefs::genHitKey((unsigned int) pad_num, 0); + if (m_deadChannelMap.contains(hitsetkey) && + std::find(m_deadChannelMap[hitsetkey].begin(), m_deadChannelMap[hitsetkey].end(), checkkey) != m_deadChannelMap[hitsetkey].end()) + { + continue; // Skip this hit, the channel is dead + } +} +if (m_maskHotChannels) +{ + TrkrDefs::hitkey checkkey = TpcDefs::genHitKey((unsigned int) pad_num, 0); + if (m_hotChannelMap.contains(hitsetkey) && + std::find(m_hotChannelMap[hitsetkey].begin(), m_hotChannelMap[hitsetkey].end(), checkkey) != m_hotChannelMap[hitsetkey].end()) + { + continue; // Skip this hit, the channel is hot + } +} + auto hitsetit = hitsetcontainer->findOrAddHitSet(hitsetkey); + auto single_hitsetit = single_hitsetcontainer->findOrAddHitSet(hitsetkey); - tpc_truth_clusterer.addhitset(hitsetkey, hitkey, neffelectrons); + for (unsigned int itb = 0; itb < adc_tbin.size(); ++itb) + { + const int tbin_num = adc_tbin[itb]; + if (tbin_num < 0 || tbin_num >= static_cast(tbins)) continue; + if (pad_num < 0 || pad_num >= static_cast(phibins)) continue; + + const double tshare = adc_tbin_share[itb]; + const float neffelectrons_bin = nelec * pad_fraction * tshare; + if (neffelectrons_bin < neffelectrons_threshold) continue; + + TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(static_cast(pad_num), static_cast(tbin_num)); + TrkrHit* hit = hitsetit->second->getHit(hitkey); + if (!hit) { hit = new TrkrHitv2(); hitsetit->second->addHitSpecificKey(hitkey, hit); } + hit->addEnergy(neffelectrons_bin); + + TrkrHit* single_hit = single_hitsetit->second->getHit(hitkey); + if (!single_hit) { single_hit = new TrkrHitv2(); single_hitsetit->second->addHitSpecificKey(hitkey, single_hit); } + single_hit->addEnergy(neffelectrons_bin); + + if (Verbosity() > 50) + { + const int layer_print = static_cast(layer_cand); + if (m_dbg_pad_hit_layer < 0 || layer_print == m_dbg_pad_hit_layer) + { + std::cout << "PadHit: layer=" << layer_cand + << " side=" << side + << " pad=" << pad_num + << " tbin=" << tbin_num + << " energy=" << neffelectrons_bin + << std::endl; + } + } - // repeat for the single_hitsetcontainer - // See if this hit already exists - TrkrHit *single_hit = nullptr; - single_hit = single_hitsetit->second->getHit(hitkey); - if (!single_hit) - { - // create a new one - single_hit = new TrkrHitv2(); - single_hitsetit->second->addHitSpecificKey(hitkey, single_hit); - } - // Either way, add the energy to it -- adc values will be added at digitization - single_hit->addEnergy(neffelectrons); + tpc_truth_clusterer.addhitset(hitsetkey, hitkey, neffelectrons_bin); - /* - if (Verbosity() > 0) - { - assert(nthit); - nthit->Fill(layernum, pad_num, tbin_num, neffelectrons); + const double tcenter = LayerGeom->get_zcenter(tbin_num); + const double phicenter = LayerGeom->get_phicenter(pad_num, side); + phi_integral += phicenter * neffelectrons_bin; + t_integral += tcenter * neffelectrons_bin; + weight += neffelectrons_bin; + } } - */ - - } // end of loop over adc T bins - } // end of loop over zigzag pads + } + } + // end of pad sharing and hit filling across candidate layers /* pass_data.phi_integral = phi_integral; */ /* pass_data.time_integral = t_integral; */ @@ -889,22 +1210,11 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( } */ - if (Verbosity() > 100) + if (Verbosity() > 100 && weight > 0) { - if (layernum == print_layer) - { - std::cout << " hit " << m_NHits << " quick centroid for this electron " << std::endl; - std::cout << " phi centroid = " << phi_integral / weight << " phi in " << phi << " phi diff " << phi_integral / weight - phi << std::endl; - std::cout << " t centroid = " << t_integral / weight << " t in " << t_gem << " t diff " << t_integral / weight - t_gem << std::endl; - // For a single track event, this captures the distribution of single electron centroids on the pad plane for layer print_layer. - // The centroid of that should match the cluster centroid found by PHG4TpcClusterizer for layer print_layer, if everything is working - // - matches to < .01 cm for a few cases that I checked - - /* - assert(nthit); - nthit->Fill(hit, layernum, phi, phi_integral / weight, t_gem, t_integral / weight, weight); - */ - } + std::cout << " hit " << m_NHits << " quick centroid for this electron " << std::endl; + std::cout << " phi centroid = " << phi_integral / weight << " phi in " << phi << " phi diff " << phi_integral / weight - phi << std::endl; + std::cout << " t centroid = " << t_integral / weight << " t in " << t_gem << " t diff " << t_integral / weight - t_gem << std::endl; } m_NHits++; @@ -929,19 +1239,19 @@ double PHG4TpcPadPlaneReadout::check_phi(const unsigned int side, const double p double daPhi = 0; if (s == 0) { - daPhi = fabs(sector_min_Phi[side][11] + 2 * M_PI - sector_max_Phi[side][s]); + daPhi = std::fabs(sector_min_Phi[side][11] + 2 * M_PI - sector_max_Phi[side][s]); } else { - daPhi = fabs(sector_min_Phi[side][s - 1] - sector_max_Phi[side][s]); + daPhi = std::fabs(sector_min_Phi[side][s - 1] - sector_max_Phi[side][s]); } double min_phi = sector_max_Phi[side][s]; double max_phi = sector_max_Phi[side][s] + daPhi; if (new_phi <= max_phi && new_phi >= min_phi) { - if (fabs(max_phi - new_phi) > fabs(new_phi - min_phi)) + if (std::fabs(max_phi - new_phi) > std::fabs(new_phi - min_phi)) { - new_phi = min_phi - phi_bin_width / 5; + new_phi = min_phi - phi_bin_width / 5; } else { @@ -958,7 +1268,6 @@ double PHG4TpcPadPlaneReadout::check_phi(const unsigned int side, const double p return new_phi; } - void PHG4TpcPadPlaneReadout::EnableSingleCloudVisualization(bool enable, const std::string &output_file, int target_side, @@ -971,7 +1280,6 @@ void PHG4TpcPadPlaneReadout::EnableSingleCloudVisualization(bool enable, m_visualization_target_layer = target_layer; m_visualization_grid_step = grid_step; m_visualization_done = false; - std::cout<<" EnableSingleCloudVisualization to "< &phibin_pad, std::vector &phibin_pad_share) +{ + const double radius = LayerGeom->get_radius(); + const double phistepsize = LayerGeom->get_phistep(); + const auto phibins = LayerGeom->get_phibins(); + + double rphi = phi * radius; + if (Verbosity() > 100) + { + if (LayerGeom->get_layer() == print_layer) + { + std::cout << " populate_zigzag_phibins for layer " << layernum << " with radius " << radius << " phi " << phi + << " rphi " << rphi << " phistepsize " << phistepsize << std::endl; + std::cout << " fcharge created: radius " << radius << " rphi " << rphi << " cloud_sig_rp " << cloud_sig_rp << std::endl; + } + } + + +} +*/ double PHG4TpcPadPlaneReadout::integratedDensityOfCircleAndPad( double hitX, double hitY, double sigma, const std::vector& pad, - double gridStep + double gridStep, + std::vector* debug_samples ) { // sanity checks if (sigma <= 0.0 || pad.empty()) return 0.0; // constants - const int nSigma = _nsigmas; + const double nSigma = _nsigmas; double R = nSigma * sigma; double gaussConst = 1.0 / (2.0 * M_PI * sigma * sigma); double expDenominator = 2.0 * sigma * sigma; @@ -1022,8 +1349,9 @@ double PHG4TpcPadPlaneReadout::integratedDensityOfCircleAndPad( } // 1) Compute pad bounding box - double minx = pad[0].x, maxx = pad.back().x; - double miny = pad[0].y, maxy = pad.back().y; + // Initialize min/max from the same first vertex to ensure both extremas consider index 0 + double minx = pad[0].x, maxx = pad[0].x; + double miny = pad[0].y, maxy = pad[0].y; for (size_t i = 1; i < pad.size(); ++i) { minx = std::min(minx, pad[i].x); maxx = std::max(maxx, pad[i].x); @@ -1032,7 +1360,7 @@ double PHG4TpcPadPlaneReadout::integratedDensityOfCircleAndPad( // std::cout<<"Pad point "<Draw("P SAME"); center_markers.push_back(center); } -std::cout<<" m_visualization_dump_file "< PHG4TpcPadPlaneReadout::getLayersToCheck( - unsigned int layernum, - double rad_gem, - double cloud_sig_rp) const -{ - std::vector layers_to_check; - - // Get current layer geometry - PHG4TpcGeom *current_layer = GeomContainer->GetLayerCellGeom(layernum); - if (!current_layer) { - layers_to_check.push_back(static_cast(layernum)); - return layers_to_check; - } - - // Calculate radial extent of the electron cloud (nsigmas coverage) - const double cloud_radial_extent = _nsigmas * cloud_sig_rp; - const double r_min = rad_gem - cloud_radial_extent; - const double r_max = rad_gem + cloud_radial_extent; - - // Get current layer boundaries - const double layer_r_center = current_layer->get_radius(); - const double layer_thickness = current_layer->get_thickness(); - const double layer_r_min = layer_r_center - layer_thickness / 2.0; - const double layer_r_max = layer_r_center + layer_thickness / 2.0; - - // Always check current layer - layers_to_check.push_back(static_cast(layernum)); - - // Determine module boundaries for valid layer range - const int tpc_module = static_cast((layernum - 7) / 16); - const int module_first_layer = 7 + tpc_module * 16; - const int module_last_layer = module_first_layer + 15; - - // Check if cloud extends into previous layer - if (r_min < layer_r_min) { - int prev_layer = static_cast(layernum) - 1; - if (prev_layer >= module_first_layer && prev_layer >= 7) { - PHG4TpcGeom *prev_geom = GeomContainer->GetLayerCellGeom(prev_layer); - if (prev_geom) { - // Verify there's actual overlap - double prev_r_max = prev_geom->get_radius() + prev_geom->get_thickness() / 2.0; - if (r_min < prev_r_max) { - layers_to_check.push_back(prev_layer); - } - } - } - } - - // Check if cloud extends into next layer - if (r_max > layer_r_max) { - int next_layer = static_cast(layernum) + 1; - if (next_layer <= module_last_layer && - next_layer < static_cast(Pads.size())) { - PHG4TpcGeom *next_geom = GeomContainer->GetLayerCellGeom(next_layer); - if (next_geom) { - // Verify there's actual overlap - double next_r_min = next_geom->get_radius() - next_geom->get_thickness() / 2.0; - if (r_max > next_r_min) { - layers_to_check.push_back(next_layer); - } - } - } - } - - return layers_to_check; -} - - -// In SERF_zigzag_phibins, add timing: -void PHG4TpcPadPlaneReadout::SERF_zigzag_phibins( - const unsigned int side, const unsigned int layernum, const double phi, - const double rad_gem, const double cloud_sig_rp, - std::vector& pad_layer, std::vector &phibin_pad, - std::vector &phibin_pad_share) +void PHG4TpcPadPlaneReadout::SERF_zigzag_phibins(const unsigned int side, const unsigned int layernum, const double phi, const double rad_gem,const double cloud_sig_rp, std::vector &phibin_pad, std::vector &phibin_pad_share) { - g_zigzag_call_count++; -/*std::cout << "[PROGRESS] SERF_zigzag_phibins call #" << g_zigzag_call_count - << " - Starting for side=" << side << ", layer=" << layernum - << ", phi=" << std::fixed << std::setprecision(4) << phi - << ", rad=" << rad_gem << std::endl;*/ - - //auto t_start = std::chrono::high_resolution_clock::now(); - + const double radius = LayerGeom->get_radius(); const double phistepsize = LayerGeom->get_phistep(); const auto phibins = LayerGeom->get_phibins(); @@ -1527,179 +1713,210 @@ void PHG4TpcPadPlaneReadout::SERF_zigzag_phibins( double xNew = 0.0; double yNew = 0.0; - rotatePointToSector(x, y, xNew, yNew, sector1); + rotatePointToSector( x, y, side, sector1, xNew, yNew); + // double phiNew = std::atan2(yNew,xNew); + //double rad_gem_new = std::sqrt(xNew*xNew+ yNew*yNew); + // std::cout<<"PHG4TpcPadPlaneReadout::SERF_zigzag_phibins: x = "<get_radius() - LayerGeom->get_thickness() / 2.0 << std::endl; + } + if (radlim_high > (LayerGeom->get_radius() + LayerGeom->get_thickness() / 2.0)) + { + std::cout << " radlim_high " << radlim_high << " is in the next layer " <get_radius() + LayerGeom->get_thickness() / 2.0 << std::endl; +} +*/ + // std::cout << " SERF zigzags: phi " << phi << " philim_low " << philim_low << " phibin_low " << phibin_low << " philim_high " << philim_high << " phibin_high " << phibin_high << " npads " << npads << std::endl; + if (capture_debug) + { + const int reserve_count = (npads >= 0) ? (npads + 2) : 2; + debug_contribs.reserve(static_cast(reserve_count)); + } - //auto t_after_phibin_calc = std::chrono::high_resolution_clock::now(); + for (int ipad = 0; ipad <= npads; ipad++) + { + int pad_now = phibin_low + ipad; + // if(phibin_low<0 && phibin_high<0) pad_now = phibin_high + ipad; + // check that we do not exceed the maximum number of pads, wrap if necessary + if (pad_now >= phibins) + { + pad_now -= phibins; + } + + int look_pad = pad_now ; + // int n = ntpc_phibins_sector[tpc_module]; + //look_pad = ( n - ( (look_pad % n) + n ) % n ) % n ; + look_pad = ((look_pad % ntpc_phibins_sector[tpc_module]) + ntpc_phibins_sector[tpc_module]) % ntpc_phibins_sector[tpc_module]; - int total_pads_processed = 0; - double total_integration_time = 0.0; - std::vector layers_to_check = getLayersToCheck(layernum, rad_gem, cloud_sig_rp); - /* + look_pad = ntpc_phibins_sector[tpc_module] - look_pad -1; - std::cout<<"LayerNum "<((check_layer - 7) / 16); - if (module_for_check != tpc_module) continue; + // std::cout<<"pad now = "< poly; + if (m_use_rectangular_pad_response) { - int pad_now = phibin_low + ipad; - if (pad_now >= phibins) pad_now -= phibins; - - int look_pad = pad_now; - while(look_pad >= ntpc_phibins_sector[tpc_module]) { - look_pad -= ntpc_phibins_sector[tpc_module]; - } - look_pad = ntpc_phibins_sector[tpc_module] - look_pad - 1; - - const auto& padinfo = Pads[check_layer][look_pad]; - - // Time the integration - auto t_int_start = std::chrono::high_resolution_clock::now(); - - double charge = integratedDensityOfCircleAndPad( - xNew, yNew, cloud_sig_rp, padinfo.vertices, 0.0, - capture_debug ? debug_samples : nullptr); - - auto t_int_end = std::chrono::high_resolution_clock::now(); - total_integration_time += std::chrono::duration(t_int_end - t_int_start).count(); - total_pads_processed++; - - if (capture_debug && charge > 0.0) + // rectangular pad approximated by radial and phi bounds + const double pad_phi_center = LayerGeom->get_phicenter(pad_now, side); + const double half_phi = 0.5 * LayerGeom->get_phistep(); + const double r_center = LayerGeom->get_radius(); + const double half_thickness = 0.5 * LayerGeom->get_thickness(); + const double r_low = r_center - half_thickness; + const double r_high = r_center + half_thickness; + + const double phi_edges[2] = {pad_phi_center - half_phi, pad_phi_center + half_phi}; + const double r_edges[2] = {r_low, r_high}; + + poly.reserve(4); + // ordered rectangle: (r_low,phi_low) -> (r_low,phi_high) -> (r_high,phi_high) -> (r_high,phi_low) + poly.push_back({r_edges[0] * std::cos(phi_edges[0]), r_edges[0] * std::sin(phi_edges[0])}); + poly.push_back({r_edges[0] * std::cos(phi_edges[1]), r_edges[0] * std::sin(phi_edges[1])}); + poly.push_back({r_edges[1] * std::cos(phi_edges[1]), r_edges[1] * std::sin(phi_edges[1])}); + poly.push_back({r_edges[1] * std::cos(phi_edges[0]), r_edges[1] * std::sin(phi_edges[0])}); + } + else + { + // Guard against missing pad polygons + if (layernum >= Pads.size() || look_pad < 0 || static_cast(look_pad) >= Pads[layernum].size()) { - DebugPadContribution dbg; - dbg.pad_bin = pad_now; - dbg.charge = charge; - dbg.pad_phi = LayerGeom->get_phicenter(pad_now, side); - dbg.polygon = padinfo.vertices; - debug_contribs.push_back(std::move(dbg)); + // No polygon data for this pad/layer; skip contribution + continue; } - - pad_layer.push_back(static_cast(check_layer)); - phibin_pad.push_back(pad_now); - phibin_pad_share.push_back(charge); + auto padinfo = Pads[layernum][look_pad]; + poly = padinfo.vertices; + } + // std::cout<<"Calculate charge for pad with cx = "<get_phicenter(pad_now, side) <<" pad phi center "<get_phicenter(pad_now, side) <<" pad phi center "<empty() && !debug_contribs.empty()) + if (capture_debug && debug_samples && !debug_contribs.empty()) { VisualizationCircle this_circle{xNew, yNew, _nsigmas * cloud_sig_rp}; - if (m_visualize_all_matches) { - m_visualization_aggregate_samples.insert( - m_visualization_aggregate_samples.end(), - debug_samples->begin(), debug_samples->end()); + m_visualization_aggregate_samples.insert(m_visualization_aggregate_samples.end(), + debug_samples->begin(), debug_samples->end()); m_visualization_circles.push_back(this_circle); - - for (const auto &pad : debug_contribs) { - m_visualization_pad_union[pad.pad_bin] = pad; - } - - maybeVisualizeAvalanche(side, layernum, phi, rad_gem, cloud_sig_rp, - xNew, yNew, debug_contribs, - m_visualization_aggregate_samples, - m_visualization_circles); + for (const auto &pad : debug_contribs) m_visualization_pad_union[pad.pad_bin] = pad; + + const double vis_x = m_use_rectangular_pad_response ? x : xNew; + const double vis_y = m_use_rectangular_pad_response ? y : yNew; + maybeVisualizeAvalanche(side, layernum, phi, rad_gem, cloud_sig_rp, vis_x, vis_y, debug_contribs, + m_visualization_aggregate_samples, m_visualization_circles); } else { std::vector single_circle{this_circle}; - maybeVisualizeAvalanche(side, layernum, phi, rad_gem, cloud_sig_rp, - xNew, yNew, debug_contribs, + const double vis_x = m_use_rectangular_pad_response ? x : xNew; + const double vis_y = m_use_rectangular_pad_response ? y : yNew; + maybeVisualizeAvalanche(side, layernum, phi, rad_gem, cloud_sig_rp, vis_x, vis_y, debug_contribs, *debug_samples, single_circle); m_visualization_done = true; } ++m_visualization_cloud_counter; } - /*auto t_end = std::chrono::high_resolution_clock::now(); - - // Print timing results - auto dur_rotation = std::chrono::duration(t_after_rotation - t_start).count(); - auto dur_viz_check = std::chrono::duration(t_after_viz_check - t_after_rotation).count(); - auto dur_alloc = std::chrono::duration(t_after_alloc - t_after_viz_check).count(); - auto dur_phibin = std::chrono::duration(t_after_phibin_calc - t_after_alloc).count(); - auto dur_total_int = std::chrono::duration(t_after_integration - t_after_phibin_calc).count(); - auto dur_viz = std::chrono::duration(t_end - t_after_integration).count(); - auto dur_total = std::chrono::duration(t_end - t_start).count(); - - std::cout << "[TIMING] SERF_zigzag_phibins side=" << side << " layer=" << layernum - << " | Total: " << dur_total << " μs" - << " | Rotation: " << dur_rotation << " μs" - << " | VizCheck: " << dur_viz_check << " μs" - << " | Alloc: " << dur_alloc << " μs" - << " | PhiBinCalc: " << dur_phibin << " μs" - << " | Integration: " << dur_total_int << " μs (" << total_pads_processed << " pads, avg " - << (total_pads_processed > 0 ? total_integration_time/total_pads_processed : 0.0) << " μs/pad)" - << " | Visualization: " << dur_viz << " μs" - << " | capture_debug=" << (capture_debug ? "true" : "false") - << std::endl; - */ + return; } -//--------------------------------------------------------- + void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int side, const unsigned int layernum, const double phi, const double cloud_sig_rp, std::vector &phibin_pad, std::vector &phibin_pad_share) { const double radius = LayerGeom->get_radius(); const double phistepsize = LayerGeom->get_phistep(); const auto phibins = LayerGeom->get_phibins(); - // make the charge distribution gaussian double rphi = phi * radius; if (Verbosity() > 100) @@ -1723,13 +1940,22 @@ void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int side, co int phibin_low = LayerGeom->get_phibin(philim_high, side); int phibin_high = LayerGeom->get_phibin(philim_low, side); int npads = phibin_high - phibin_low; + int sector = 0; + for (int i=0;i<12;i++) + { + if (phi >= sector_min_Phi[side][i] && phi < sector_max_Phi[side][i]) + { + sector = i; + break; + } + } if (Verbosity() > 1000) { if (layernum == print_layer) { - std::cout << " zigzags: phi " << phi << " philim_low " << philim_low << " phibin_low " << phibin_low - << " philim_high " << philim_high << " phibin_high " << phibin_high << " npads " << npads << std::endl; + std::cout << " REF zigzags: phi " << phi << " philim_low " << philim_low << " phibin_low " << phibin_low + << " philim_high " << philim_high << " phibin_high " << phibin_high << " npads " << npads <<" sector "<get_phicenter(pad_now, side); sum_of_pads_phi += pads_phi[ipad]; - sum_of_pads_absphi += fabs(pads_phi[ipad]); + sum_of_pads_absphi += std::fabs(pads_phi[ipad]); } for (int ipad = 0; ipad <= npads; ipad++) @@ -1800,7 +2026,7 @@ void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int side, co const double sigma = cloud_sig_rp; // eshulga // Checking if the pads are on the same side of the TPC in phi - if (fabs(sum_of_pads_phi) != sum_of_pads_absphi) + if (std::fabs(sum_of_pads_phi) != sum_of_pads_absphi) { if (phi < -M_PI / 2 && phi_now > 0) { @@ -1812,7 +2038,7 @@ void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int side, co } if (phi < 0 && phi_now > 0) { - x_loc_tmp = (phi_now + fabs(phi)) * radius; + x_loc_tmp = (phi_now + std::fabs(phi)) * radius; } if (phi > 0 && phi_now < 0) { @@ -1821,13 +2047,24 @@ void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int side, co } const double x_loc = x_loc_tmp; - // calculate fraction of the total charge on this strip - /* - this corresponds to integrating the charge distribution Gaussian function (centered on rphi and of width cloud_sig_rp), - convoluted with a strip response function, which is triangular from -pitch to +pitch, with a maximum of 1. at stript center - */ - overlap[ipad] = - (pitch - x_loc) * (std::erf(x_loc / (M_SQRT2 * sigma)) - std::erf((x_loc - pitch) / (M_SQRT2 * sigma))) / (pitch * 2) + (pitch + x_loc) * (std::erf((x_loc + pitch) / (M_SQRT2 * sigma)) - std::erf(x_loc / (M_SQRT2 * sigma))) / (pitch * 2) + (gaus(x_loc - pitch, sigma) - gaus(x_loc, sigma)) * square(sigma) / pitch + (gaus(x_loc + pitch, sigma) - gaus(x_loc, sigma)) * square(sigma) / pitch; + // calculate fraction of the total charge on this pad + if (m_use_rectangular_pad_response) + { + // flat response across the pad width (rectangular pad) + const double half_width = pad_rphi / 2.0; + const double upper = (x_loc + half_width) / (M_SQRT2 * sigma); + const double lower = (x_loc - half_width) / (M_SQRT2 * sigma); + overlap[ipad] = 0.5 * (std::erf(upper) - std::erf(lower)); + } + else + { + /* + this corresponds to integrating the charge distribution Gaussian function (centered on rphi and of width cloud_sig_rp), + convoluted with a strip response function, which is triangular from -pitch to +pitch, with a maximum of 1. at strip center + */ + overlap[ipad] = + (pitch - x_loc) * (std::erf(x_loc / (M_SQRT2 * sigma)) - std::erf((x_loc - pitch) / (M_SQRT2 * sigma))) / (pitch * 2) + (pitch + x_loc) * (std::erf((x_loc + pitch) / (M_SQRT2 * sigma)) - std::erf(x_loc / (M_SQRT2 * sigma))) / (pitch * 2) + (gaus(x_loc - pitch, sigma) - gaus(x_loc, sigma)) * square(sigma) / pitch + (gaus(x_loc + pitch, sigma) - gaus(x_loc, sigma)) * square(sigma) / pitch; + } } // now we have the overlap for each pad @@ -1835,11 +2072,15 @@ void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int side, co { phibin_pad.push_back(pad_keep[ipad]); phibin_pad_share.push_back(overlap[ipad]); + std::cout<<" REF zigzags: ipad " << ipad << " pad_now " << pad_keep[ipad] << " charge " << overlap[ipad] + << " phi center " << LayerGeom->get_phicenter(pad_keep[ipad], side) << std::endl; } return; } + + void PHG4TpcPadPlaneReadout::UseGain(const int flagToUseGain) { m_flagToUseGain = flagToUseGain; @@ -1876,8 +2117,13 @@ void PHG4TpcPadPlaneReadout::SetDefaultParameters() set_default_double_param("tpc_maxradius_inner", 40.249); // 40.0); // cm set_default_double_param("tpc_maxradius_mid", 57.475); // 60.0); set_default_double_param("tpc_maxradius_outer", 75.911); // 77.0); // from Tom + set_default_double_param("tpc_sampa_peaking_time", 80.0); // ns - set_default_double_param("neffelectrons_threshold", 1.0); + // Minimum effective electrons per (pad,tbin) to create a hit + // Set to 0.0 by default to include all contributions + set_default_double_param("neffelectrons_threshold", 0.0); + set_default_double_param("maxdriftlength", 105.5); // cm + set_default_double_param("tpc_adc_clock", 53.326184); // ns, for 18.8 MHz clock set_default_double_param("gem_cloud_sigma", 0.04); // cm = 400 microns set_default_double_param("sampa_shaping_lead", 32.0); // ns, for 80 ns SAMPA set_default_double_param("sampa_shaping_tail", 48.0); // ns, for 80 ns SAMPA @@ -1918,11 +2164,20 @@ void PHG4TpcPadPlaneReadout::UpdateInternalParameters() sigmaL = {{get_double_param("sampa_shaping_lead"), get_double_param("sampa_shaping_tail")}}; + const double tpc_adc_clock = get_double_param("tpc_adc_clock"); + + const double MaxZ = get_double_param("maxdriftlength"); + const double TBinWidth = tpc_adc_clock; + const double MaxT = extended_readout_time + 2.0 * MaxZ / drift_velocity; // allows for extended time readout + const double MinT = 0; + NTBins = (int) ((MaxT - MinT) / TBinWidth) + 1; + averageGEMGain = get_double_param("gem_amplification"); polyaTheta = get_double_param("polya_theta"); -} + Ts = get_double_param("tpc_sampa_peaking_time"); +} void PHG4TpcPadPlaneReadout::makeChannelMask(hitMaskTpc &aMask, const std::string &dbName, const std::string &totalChannelsToMask) { CDBTTree *cdbttree; @@ -1958,9 +2213,13 @@ void PHG4TpcPadPlaneReadout::makeChannelMask(hitMaskTpc &aMask, const std::strin } delete cdbttree; + } +// ------------------------------------------------------------------------- +// REPLACEMENT FUNCTIONS (From Code A) +// ------------------------------------------------------------------------- -void PHG4TpcPadPlaneReadout::sampaTimeDistribution(double tzero, std::vector &adc_tbin, std::vector &adc_tbin_share) +void PHG4TpcPadPlaneReadout::sampaTimeDistribution(double tzero, std::vector &adc_tbin, std::vector &adc_tbin_share) { // tzero is the arrival time of the electron at the GEM // Ts is the sampa peaking time @@ -1972,64 +2231,48 @@ void PHG4TpcPadPlaneReadout::sampaTimeDistribution(double tzero, std::vectorget_zcenter(tbinzero) + tstepsize/2.0; - double vfirst_end = sampaShapingResponseFunction(tzero, tfirst_end); + double vfirst_end = sampaShapingResponseFunction(tzero, tfirst_end); double first_integral = (vfirst_end / 2.0) * (tfirst_end - tzero); adc_tbin.push_back(tbinzero); adc_tbin_share.push_back(first_integral); - /* - if (LayerGeom->get_layer() == print_layer) - { - std::cout << " tzero " << tzero << " tbinzero " << tbinzero << " iclock 0 " - << " tfirst_end " << tfirst_end << " vfirst_end " << vfirst_end << " first_integral " << first_integral << std::endl; - } - */ - for(int iclock = 1; iclock < nclocks; ++iclock) + { + int tbin = tbinzero + iclock; + if (tbin < 0 || tbin > LayerGeom->get_zbins()) { - int tbin = tbinzero + iclock; - if (tbin < 0 || tbin > LayerGeom->get_zbins()) - { - if (Verbosity() > 0) + if (Verbosity() > 0) { std::cout << " t bin " << tbin << " is outside range of " << LayerGeom->get_zbins() << " so skip it" << std::endl; } - continue; - } - - // get the beginning and end of this clock bin - double tcenter = LayerGeom->get_zcenter(tbin); - double tlow = tcenter - tstepsize/2.0; - - // sample the voltage in this bin at nsamples-1 locations - int nsamples = 6; - double sample_step = tstepsize / (double) nsamples; - double sintegral = 0; - for(int isample = 0; isample < nsamples; ++isample) - { - double tnow = tlow + (double) isample * sample_step + sample_step / 2.0; - double vnow = sampaShapingResponseFunction(tzero, tnow); - sintegral += vnow * sample_step; + continue; + } - /* - if (LayerGeom->get_layer() == print_layer) - { - std::cout << " tzero " << tzero << " tbinzero " << tbinzero << " iclock " << iclock << " tbin " << tbin << " isample " << isample - << " tnow " << tnow << " vnow " << vnow << " sintegral " << sintegral << std::endl; - } - */ - } + // get the beginning and end of this clock bin + double tcenter = LayerGeom->get_zcenter(tbin); + double tlow = tcenter - tstepsize/2.0; - - adc_tbin.push_back(tbin); - adc_tbin_share.push_back(sintegral); + // sample the voltage in this bin at nsamples locations + int nsamples = 6; + double sample_step = tstepsize / (double) nsamples; + double sintegral = 0; + for(int isample = 0; isample < nsamples; ++isample) + { + double tnow = tlow + (double) isample * sample_step + sample_step / 2.0; + double vnow = sampaShapingResponseFunction(tzero, tnow); + sintegral += vnow * sample_step; } + + adc_tbin.push_back(tbin); + adc_tbin_share.push_back(sintegral); + } } double PHG4TpcPadPlaneReadout::sampaShapingResponseFunction(double tzero, double t) const - { - double v = exp(-4*(t-tzero)/Ts) * pow( (t-tzero)/Ts, 4.0); - - return v; - } +{ + // The specific SAMPA response function: V ~ (t/tau)^4 * exp(-4t/tau) + //if (t < tzero) return 0.0; + double v = exp(-4.0 * (t - tzero) / Ts) * pow((t - tzero) / Ts, 4.0); + return v; +} diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h index d083447be6..afe64de980 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h @@ -14,7 +14,6 @@ #include // for string #include #include - typedef std::map> hitMaskTpc; class PHCompositeNode; @@ -36,22 +35,32 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane int InitRun(PHCompositeNode *topNode) override; void UseGain(const int flagToUseGain); - void SetUseModuleGainWeights(const int flag) { m_use_module_gain_weights = flag; } - void SetModuleGainWeightsFileName(const std::string &name) { m_tpc_module_gain_weights_file = name; } + void SetUseModuleGainWeights(const int flag) {m_use_module_gain_weights = flag;} + void SetModuleGainWeightsFileName(const std::string &name) {m_tpc_module_gain_weights_file = name;} void ReadGain(); - void SetUsePolyaGEMGain(const int flagPolya) { m_usePolya = flagPolya; } - void SetUseLangauGEMGain(const int flagLangau) { m_useLangau = flagLangau; } - void SetLangauParsFileName(const std::string &name) { m_tpc_langau_pars_file = name; } - + void SetUsePolyaGEMGain(const int flagPolya) {m_usePolya = flagPolya;} + void SetUseLangauGEMGain(const int flagLangau) {m_useLangau = flagLangau;} + void SetLangauParsFileName(const std::string &name) {m_tpc_langau_pars_file = name;} + // Pad-sharing method selection: true → SERF polygon overlap; false → analytic triangle + void UseSerfPadSharing(bool use_serf) { m_use_serf_padsharing = use_serf; } + // If true and SERF polygons are unavailable, abort InitRun with an error + void RequireSerfPadSharing(bool require) { m_require_serf = require; } + // If true, also require polygons for all readout layers (partial coverage aborts) + void RequireSerfFullCoverage(bool require) { m_require_serf_full = require; } + // When using analytic sharing (no SERF polygons), choose rectangular response instead of triangular + void UseRectangularPadResponse(bool use_rectangular) { m_use_rectangular_pad_response = use_rectangular; } + + void SetDriftVelocity(double vd) override { drift_velocity = vd; } + void SetReadoutTime(float t) override { extended_readout_time = t; } // otherwise warning of inconsistent overload since only one MapToPadPlane methow is overridden using PHG4TpcPadPlane::MapToPadPlane; - void MapToPadPlane(TpcClusterBuilder &tpc_truth_clusterer, TrkrHitSetContainer *single_hitsetcontainer, TrkrHitSetContainer *hitsetcontainer, TrkrHitTruthAssoc * /*hittruthassoc*/, const double x_gem, const double y_gem, const double t_gem, const unsigned int side, PHG4HitContainer::ConstIterator hiter, TNtuple * /*ntpad*/, TNtuple * /*nthit*/) override; + void MapToPadPlane(TpcClusterBuilder &tpc_clustbuilder, TrkrHitSetContainer *single_hitsetcontainer, TrkrHitSetContainer *hitsetcontainer, TrkrHitTruthAssoc * /*hittruthassoc*/, const double x_gem, const double y_gem, const double t_gem, const unsigned int side, PHG4HitContainer::ConstIterator hiter, TNtuple * /*ntpad*/, TNtuple * /*nthit*/ ) override; + //void MapToPadPlane(TpcClusterBuilder &tpc_clustbuilder, TrkrHitSetContainer *single_hitsetcontainer, TrkrHitSetContainer *hitsetcontainer, TrkrHitTruthAssoc * /*hittruthassoc*/, const double x_gem, const double y_gem, const double t_gem, const unsigned int side, PHG4HitContainer::ConstIterator hiter, TNtuple * /*ntpad*/, TNtuple * /*nthit*/, TH2* h_adc_ref, TH2* h_adc_serf ) override; void SetDefaultParameters() override; void UpdateInternalParameters() override; - - void SetMaskChannelsFromFile() +void SetMaskChannelsFromFile() { m_maskFromFile = true; } @@ -65,10 +74,9 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane m_maskHotChannels = true; m_hotChannelMapName = hmap; } - void LoadAllPadPlanes(); - // Debug printing helpers + // Debug printing helpers // If set >= 0, limit PadHit prints to a single layer number; otherwise prints for all layers. - // void SetDebugPadHitLayer(int layer) { m_dbg_pad_hit_layer = layer; } + void SetDebugPadHitLayer(int layer) { m_dbg_pad_hit_layer = layer; } // Enable a one-shot visualization of a single avalanche cloud overlap with zigzag pads. // Passing target_side/target_layer < 0 matches the first cloud encountered. // grid_step <= 0 defaults to sigma/30 sampling. @@ -80,51 +88,67 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane void SetVisualizationDumpFile(const std::string &file); void SetVisualizeAllClouds(bool enable); + protected: + double Ts = 80.0; // SAMPA peaking time + + double sampaShapingResponseFunction(double tzero, double t) const; + + void sampaTimeDistribution(double tzero, + std::vector &adc_tbin, + std::vector &adc_tbin_share); + private: + // void populate_rectangular_phibins(const unsigned int layernum, const double phi, const double cloud_sig_rp, std::vector &pad_phibin, std::vector &pad_phibin_share); - void SERF_zigzag_phibins(const unsigned int side, const unsigned int layernum, const double phi, const double rad_gem, const double cloud_sig_rp, std::vector& pad_layer, std::vector &pad_phibin, std::vector &pad_phibin_share); - void populate_zigzag_phibins(const unsigned int side, const unsigned int layernum, const double phi, const double cloud_sig_rp, std::vector &phibin_pad, std::vector &phibin_pad_share); + void populate_zigzag_phibins(const unsigned int side, const unsigned int layernum, const double phi, const double cloud_sig_rp, std::vector &pad_phibin, std::vector &pad_phibin_share); + void SERF_zigzag_phibins(const unsigned int side, const unsigned int layernum, const double phi, const double rad_gem, const double cloud_sig_rp, std::vector &pad_phibin, std::vector &pad_phibin_share); + //void populate_tbins(const double t, const std::array &cloud_sig_tt, std::vector &adc_tbin, std::vector &adc_tbin_share); - void sampaTimeDistribution(double tzero, std::vector &adc_tbin, std::vector &adc_tbin_share); - double sampaShapingResponseFunction(double tzero, double t) const; - double check_phi(const unsigned int side, const double phi, const double radius); +void makeChannelMask(hitMaskTpc& aMask, const std::string& dbName, const std::string& totalChannelsToMask); + // utility: pick layers whose annulus intersects a radial window around rad + std::vector layersInRadialWindow(double rad, double sigma, double nsig) const; + + // utility: find geometry for a given layer + PHG4TpcGeom* getGeomForLayer(unsigned int layer) const; - void makeChannelMask(hitMaskTpc& aMask, const std::string& dbName, const std::string& totalChannelsToMask); + // utility: determine sector for (x,y) and rotate to a canonical frame + void rotatePointToSector(double x, double y, unsigned int side, int& sectorFound, double& xNew, double& yNew); PHG4TpcGeomContainer *GeomContainer = nullptr; PHG4TpcGeom *LayerGeom = nullptr; - double neffelectrons_threshold {std::numeric_limits::quiet_NaN()}; + // Minimum effective electrons per (pad,tbin) to create a hit. + // Default to 0.0 so all contributions are kept unless configured otherwise. + double neffelectrons_threshold = 0.0; std::array MinRadius{}; std::array MaxRadius{}; - static constexpr int NSides {2}; - static constexpr int NSectors {12}; - static const int NRSectors {3}; + static constexpr int NSides = 2; + static constexpr int NSectors = 12; + static const int NRSectors = 3; - double sigmaT {std::numeric_limits::quiet_NaN()}; + double sigmaT = std::numeric_limits::signaling_NaN(); std::array sigmaL{}; double phi_bin_width{}; - - int NTBins {std::numeric_limits::max()}; - int m_NHits {0}; + double drift_velocity = 8.0e-03; // default value, override from macro + float extended_readout_time = 0; // ns + int NTBins = std::numeric_limits::max(); + int m_NHits = 0; // Using Gain maps is turned off by default - int m_flagToUseGain {0}; + int m_flagToUseGain = 0; // Optionally apply a module-by-module weight to the GEM gain // Weights are input from a file for all 72 TPC modules - bool m_use_module_gain_weights {false}; - std::string m_tpc_module_gain_weights_file; + bool m_use_module_gain_weights = false; + std::string m_tpc_module_gain_weights_file = ""; // gaussian sampling - static constexpr double _nsigmas {5}; - - double Ts {55.0}; // SAMPA v5 peaking time + static constexpr double _nsigmas = 2.5; - double averageGEMGain {std::numeric_limits::quiet_NaN()}; - double polyaTheta {std::numeric_limits::quiet_NaN()}; + double averageGEMGain = std::numeric_limits::signaling_NaN(); + double polyaTheta = std::numeric_limits::signaling_NaN(); std::array, NSides> sector_min_Phi; std::array, NSides> sector_max_Phi; @@ -132,35 +156,33 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane // return random distribution of number of electrons after amplification of GEM for each initial ionizing electron double getSingleEGEMAmplification(); double getSingleEGEMAmplification(double weight); - static double getSingleEGEMAmplification(TF1 *f); - bool m_usePolya {false}; + double getSingleEGEMAmplification(TF1 *f); + bool m_usePolya = false; - bool m_useLangau {false}; - std::string m_tpc_langau_pars_file; + bool m_useLangau = false; + std::string m_tpc_langau_pars_file = ""; - gsl_rng *RandomGenerator {nullptr}; + gsl_rng *RandomGenerator = nullptr; std::array h_gain{nullptr}; - double m_module_gain_weight[2][3][12] { - {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, - {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, - {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}}, - {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, - {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, - {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}}}; - - TF1 *flangau[2][3][12] {{{nullptr}}}; + double m_module_gain_weight[2][3][12] = { + { {1,1,1,1,1,1,1,1,1,1,1,1}, + {1,1,1,1,1,1,1,1,1,1,1,1}, + {1,1,1,1,1,1,1,1,1,1,1,1} }, + { {1,1,1,1,1,1,1,1,1,1,1,1}, + {1,1,1,1,1,1,1,1,1,1,1,1}, + {1,1,1,1,1,1,1,1,1,1,1,1} } + }; - hitMaskTpc m_deadChannelMap; + TF1 *flangau[2][3][12] = {}; + hitMaskTpc m_deadChannelMap; hitMaskTpc m_hotChannelMap; - - bool m_maskDeadChannels {false}; +bool m_maskDeadChannels {false}; bool m_maskHotChannels {false}; bool m_maskFromFile {false}; std::string m_deadChannelMapName; std::string m_hotChannelMapName; - void loadPadPlanes(); struct Point { double x, y; }; struct PadInfo { @@ -178,31 +200,25 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane cx = cy = rad = phi = 0.0; vertices.clear(); isedge = false; - } + } }; - - struct DebugSample + + struct DebugSample { double x = 0.0; double y = 0.0; double density = 0.0; }; - - -int findPadForPoint( double x, double y, int tpc_module); +std::array,3*16+7> Pads; bool pointInPolygon( double x, double y,const std::vector& poly); -//double integratedDensityOfCircleAndPad(double hitX,double hitY, double sigma, const std::vector& pad,double gridStep = 0.0); double integratedDensityOfCircleAndPad(double hitX,double hitY, double sigma, const std::vector& pad,double gridStep = 0.0, std::vector* debug_samples = nullptr); - -const std::vector brdMaps_ = { - "/sphenix/user/mitrankova/Simulation/PadPlane/AutoPad-R1-RevA.brd", - "/sphenix/user/mitrankova/Simulation/PadPlane/AutoPad-R2-RevA-Pads.brd", - "/sphenix/user/mitrankova/Simulation/PadPlane/AutoPad-R3-RevA.brd" -}; -std::array,3*16+7> Pads; + // hard‑coded list of input .brd files + static const std::vector brdMaps_; +void loadPadPlanes(); int ntpc_phibins_sector[3] = { 94, 128, 192 }; +int findPadForPoint( double x, double y, int tpc_module); const std::array Thickness = {{ 0.56598621677629212, @@ -214,12 +230,25 @@ int ntpc_phibins_sector[3] = { 94, 128, 192 }; double min_radii_module[3]={314.9836110818037, 416.59202613529567, 589.1096495597712}; double max_radii_module[3]={399.85222874031024, 569.695373910603, 753.6667758418596}; + // choose between SERF polygon overlap (default) and triangle response + bool m_use_serf_padsharing = true; + bool m_require_serf = false; + bool m_require_serf_full = false; + bool m_serf_polygons_present = false; // summarized availability after InitRun + bool m_warned_serf_fallback = false; // printed once if per-hit fallback occurs + // analytic sharing option: rectangular (flat) pad response instead of triangular + bool m_use_rectangular_pad_response = false; + + // Debug: restrict PadHit prints to a single layer (or all if < 0) + int m_dbg_pad_hit_layer = 49; + struct DebugPadContribution { int pad_bin = -1; double charge = 0.0; std::vector polygon; double pad_phi = 0.0; + double pad_r = 0.0; }; struct VisualizationCircle { @@ -227,9 +256,7 @@ double max_radii_module[3]={399.85222874031024, 569.695373910603, 753.6667758418 double y = 0.0; double radius = 0.0; }; - size_t g_map_call_count = 0; - size_t g_zigzag_call_count = 0; - size_t g_integration_call_count = 0; + void maybeVisualizeAvalanche(unsigned int side, unsigned int layernum, double phi, @@ -240,7 +267,6 @@ double max_radii_module[3]={399.85222874031024, 569.695373910603, 753.6667758418 const std::vector &contribs, const std::vector &samples, const std::vector &circles); - std::vector getLayersToCheck(unsigned int layernum, double rad_gem, double cloud_sig_rp) const; bool m_visualize_single_cloud = false; bool m_visualization_done = false; @@ -256,6 +282,7 @@ double max_radii_module[3]={399.85222874031024, 569.695373910603, 753.6667758418 std::vector m_visualization_circles; std::map m_visualization_pad_union; + }; #endif