From 7e0fa17644ea3457e3eff18601eadb3501754d5b Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Wed, 23 Jul 2025 11:20:12 -0400 Subject: [PATCH 01/13] SERF first attempts. Brd is read --- .../g4tpc/PHG4TpcPadPlaneReadout.cc | 261 +++++++++++++++++- .../g4tpc/PHG4TpcPadPlaneReadout.h | 39 ++- 2 files changed, 297 insertions(+), 3 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index 867088964c..9041cad754 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -190,12 +190,170 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) } } } + std::cout<<"!!!!!Before Load Maps"< + 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::LoadAllPadPlanes() +{ + std::cout<<"!!!!!LoadAllPadPlanes"<& allVertices, + std::vector& allNames +) { + std::ifstream in(filename.c_str()); + if (!in) { + std::cerr << "Cannot open " << filename << "\n"; + return; + } +std::cout<<"!!!!!getPadCoordinates filename "< but haven't hit its yet? + if ( inSignal && keepSignal && !inPolygon && line.find("") == 0) { + allVertices.push_back(currentVerts); + allNames .push_back(currentName); + inPolygon = false; + continue; + } + + // 5) end of this signal + if (inSignal && line.find("") == 0) { + inSignal = false; + keepSignal = false; + continue; + } + + // 6) anything else in the file is ignored + } +} + + +std::vector +PHG4TpcPadPlaneReadout::processPadVertices( + const std::vector& vertices, + const std::vector& names +) { + std::vector centroids; + centroids.reserve(vertices.size()); + + for (size_t i = 0; i < vertices.size(); ++i) { + double sumX = 0, sumY = 0; + for (size_t j = 0; j < vertices[i].size(); ++j) { + sumX += vertices[i][j].x; + sumY += vertices[i][j].y; + } + double cx = sumX / vertices[i].size(); + double cy = sumY / vertices[i].size(); + centroids.push_back(PadCentroid{ names[i], cx, cy }); + } + return centroids; +} + +std::vector> +PHG4TpcPadPlaneReadout::loadPadPlanes( + const std::vector& filenames +) { + std::cout<<"!!!!!loadPadPlanes"<> allModules; + allModules.reserve(filenames.size()); + std::cout<<"!!!!!filenames.size() "< verts; + std::vector names; + getPadCoordinates(filenames[i], verts, names); + + std::vector module = + processPadVertices(verts, names); + + std::cout << "Module " << i + << ": Loaded " << module.size() << " pads\n"; + + allModules.push_back(module); + } + return allModules; +} + //_________________________________________________________ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification() { @@ -271,6 +429,62 @@ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification(TF1 *f) } +//_________________________________________________________ +inline void rotatePointToSector(double x, double y, + double& xNew, double& yNew, + int& sector) +{ + // 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); +} +//_________________________________________________________ +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))); +} + +//_________________________________________________________ + + + + + void PHG4TpcPadPlaneReadout::MapToPadPlane( TpcClusterBuilder &tpc_truth_clusterer, @@ -708,6 +922,51 @@ double PHG4TpcPadPlaneReadout::check_phi(const unsigned int side, const double p return new_phi; } +void PHG4TpcPadPlaneReadout::build_serf_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(); + + 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; + } + } + + +} + +void getZigzagPadFractions( + double x0_mm, + double y0_mm, + const unsigned int layer, + double sigma_r, + std::map pad_fractions + ) { + double r0 = std::hypot(x0_mm, y0_mm); + int module = (layer - 7)/16; + if(module < 0) return {{}, module}; + + auto hit_pads = edgePadsHit(x0_mm, y0_mm, sigma_r, df, 100); + + for(int pad_id : hit_pads) { + const PadRow& row = *std::find_if(df.begin(), df.end(), + [&](const PadRow& pr){ return pr.PadNumber==pad_id; }); + double fraction = integratedDensityOfCircleAndPad( + x0_mm, y0_mm, sigma_r, row.PadPath); + if(fraction > 1e-8) pad_fractions[pad_id] = fraction; + } + //double total = 0; + //for(auto& kv : pad_fractions) total += kv.second; + //if(total>0) for(auto& kv : pad_fractions) kv.second /= total; + +} 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) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h index 1c98657043..c80a5b8525 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h @@ -49,8 +49,14 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane void SetDefaultParameters() override; void UpdateInternalParameters() override; - + void LoadAllPadPlanes(); + private: + struct PadCentroid; + public: + const std::vector>& GetCentroids() const + { return centroids_; } 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 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 populate_tbins(const double t, const std::array &cloud_sig_tt, std::vector &adc_tbin, std::vector &adc_tbin_share); @@ -117,7 +123,36 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane TF1 *flangau[2][3][12] = {{{nullptr}}}; - + struct Point { double x, y; }; + + // one pad’s name + centroid + struct PadCentroid { + std::string name; + double cx, cy; + }; + + // a pad’s polygon + typedef std::vector PadVertices; + + // hard‑coded list of input .brd files + static const std::vector brdMaps_; + + // parse one .brd into vertices+names + void getPadCoordinates(const std::string& filename, + std::vector& allVertices, + std::vector& allNames); + + // compute centroids from those raw vertices + std::vector + processPadVertices(const std::vector& vertices, + const std::vector& names); + + // drive everything: loop over brdMaps_, call getPadCoordinates+process… + std::vector> + loadPadPlanes(const std::vector& filenames); + + // storage for the results + std::vector> centroids_; }; #endif From 682ce0dd92ed852f1631cccef000d4e340431c3a Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Thu, 24 Jul 2025 15:32:53 -0400 Subject: [PATCH 02/13] Correct read of .brd to [mod][pad]{properties} --- .../g4tpc/PHG4TpcPadPlaneReadout.cc | 144 ++++++------------ .../g4tpc/PHG4TpcPadPlaneReadout.h | 40 ++--- 2 files changed, 61 insertions(+), 123 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index 9041cad754..c8624c5b9b 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -191,17 +191,25 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) } } std::cout<<"!!!!!Before Load Maps"< "/sphenix/user/mitrankova/Simulation/PadPlane/AutoPad-R3-RevA.brd" }; -void PHG4TpcPadPlaneReadout::LoadAllPadPlanes() -{ - std::cout<<"!!!!!LoadAllPadPlanes"<& allVertices, - std::vector& allNames -) { - std::ifstream in(filename.c_str()); +void PHG4TpcPadPlaneReadout::loadPadPlanes() { + for (size_t i = 0; i < brdMaps_.size(); ++i) { + + std::ifstream in(brdMaps_[i].c_str()); if (!in) { - std::cerr << "Cannot open " << filename << "\n"; + std::cerr << "Cannot open " << brdMaps_[i] << "\n"; return; } -std::cout<<"!!!!!getPadCoordinates filename "< but haven't hit its yet? if ( inSignal && keepSignal && !inPolygon && line.find("") == 0) { - allVertices.push_back(currentVerts); - allNames .push_back(currentName); inPolygon = false; + p.cx = sumX / p.vertices.size(); + p.cy = sumY / p.vertices.size(); + p.id=iter; + iter++; continue; } - // 5) end of this signal + if (inSignal && line.find("") == 0) { inSignal = false; keepSignal = false; + + Pads[i].push_back(p); + continue; } - // 6) anything else in the file is ignored } } - - -std::vector -PHG4TpcPadPlaneReadout::processPadVertices( - const std::vector& vertices, - const std::vector& names -) { - std::vector centroids; - centroids.reserve(vertices.size()); - - for (size_t i = 0; i < vertices.size(); ++i) { - double sumX = 0, sumY = 0; - for (size_t j = 0; j < vertices[i].size(); ++j) { - sumX += vertices[i][j].x; - sumY += vertices[i][j].y; - } - double cx = sumX / vertices[i].size(); - double cy = sumY / vertices[i].size(); - centroids.push_back(PadCentroid{ names[i], cx, cy }); - } - return centroids; -} - -std::vector> -PHG4TpcPadPlaneReadout::loadPadPlanes( - const std::vector& filenames -) { - std::cout<<"!!!!!loadPadPlanes"<> allModules; - allModules.reserve(filenames.size()); - std::cout<<"!!!!!filenames.size() "< verts; - std::vector names; - getPadCoordinates(filenames[i], verts, names); - - std::vector module = - processPadVertices(verts, names); - - std::cout << "Module " << i - << ": Loaded " << module.size() << " pads\n"; - - allModules.push_back(module); - } - return allModules; } //_________________________________________________________ @@ -428,7 +382,7 @@ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification(TF1 *f) return nelec; } - +/* //_________________________________________________________ inline void rotatePointToSector(double x, double y, double& xNew, double& yNew, @@ -481,7 +435,7 @@ double gaussianIntegral1D(double a, double b, double mu, double sigma) { //_________________________________________________________ - +*/ @@ -922,6 +876,7 @@ double PHG4TpcPadPlaneReadout::check_phi(const unsigned int side, const double p return new_phi; } +/* void PHG4TpcPadPlaneReadout::build_serf_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(); @@ -941,7 +896,8 @@ void PHG4TpcPadPlaneReadout::build_serf_zigzag_phibins(const unsigned int side, } - +*/ +/* void getZigzagPadFractions( double x0_mm, double y0_mm, @@ -966,7 +922,7 @@ void getZigzagPadFractions( //for(auto& kv : pad_fractions) total += kv.second; //if(total>0) for(auto& kv : pad_fractions) kv.second /= total; -} +}*/ 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) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h index c80a5b8525..54aadd304a 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h @@ -49,16 +49,13 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane void SetDefaultParameters() override; void UpdateInternalParameters() override; - void LoadAllPadPlanes(); - private: - struct PadCentroid; - public: - const std::vector>& GetCentroids() const - { return centroids_; } + + 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 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 build_serf_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 populate_tbins(const double t, const std::array &cloud_sig_tt, std::vector &adc_tbin, std::vector &adc_tbin_share); double check_phi(const unsigned int side, const double phi, const double radius); @@ -125,34 +122,19 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane struct Point { double x, y; }; - // one pad’s name + centroid - struct PadCentroid { - std::string name; - double cx, cy; + struct PadInfo { + std::string name; // pad name + int id; // pad id + double cx, cy; // centroid coords + std::vector vertices; // pad polygon }; - - // a pad’s polygon - typedef std::vector PadVertices; + +std::array,3> Pads; // hard‑coded list of input .brd files static const std::vector brdMaps_; +void loadPadPlanes(); - // parse one .brd into vertices+names - void getPadCoordinates(const std::string& filename, - std::vector& allVertices, - std::vector& allNames); - - // compute centroids from those raw vertices - std::vector - processPadVertices(const std::vector& vertices, - const std::vector& names); - - // drive everything: loop over brdMaps_, call getPadCoordinates+process… - std::vector> - loadPadPlanes(const std::vector& filenames); - - // storage for the results - std::vector> centroids_; }; #endif From bfc458bcc91a50f789346fc1a2c92559653a8770 Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Fri, 1 Aug 2025 17:13:59 -0400 Subject: [PATCH 03/13] Read Pads and calculate integral. Don't know if works --- .../g4tpc/PHG4TpcPadPlaneReadout.cc | 291 +++++++++++++++--- .../g4tpc/PHG4TpcPadPlaneReadout.h | 38 ++- 2 files changed, 277 insertions(+), 52 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index c8624c5b9b..8161e5f191 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -201,7 +201,7 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) std::cout << "[PHG4TpcPadPlaneReadout] Module " << i << " has " << Pads[i].size() << " pads\n"; } - int check_pad = 0; + int check_pad = 81; std::cout<<"Pad in module 0, pad 0: " << Pads[0][check_pad].name << " at (" << Pads[0][check_pad].cx << ", " << Pads[0][check_pad].cy << ")\n"; @@ -237,10 +237,16 @@ void PHG4TpcPadPlaneReadout::loadPadPlanes() { bool inSignal = false; bool inPolygon = false; bool keepSignal = false; - PadInfo p; + std::string pname_tmp; int iter = 0; double sumX = 0, sumY = 0; + PadInfo p; + int layer = -1, pad = -1; + p.clear(); + p.isedge = false; + for(int l=0;l<7;l++) + Pads[l].push_back(p); while (std::getline(in, line)) { @@ -264,7 +270,17 @@ void PHG4TpcPadPlaneReadout::loadPadPlanes() { } if ( inSignal && keepSignal && !inPolygon && line.find("") == 0) { inSignal = false; keepSignal = false; - - Pads[i].push_back(p); + std::cout<<"Module "< computeShaperKernel() { int NT = static_cast(DetectorParams::window_ns / DetectorParams::adc_dt); @@ -436,9 +462,50 @@ double gaussianIntegral1D(double a, double b, double mu, double sigma) { //_________________________________________________________ */ - - - +/* +bool PHG4TpcPadPlaneReadout::pointInPolygon( + double x, double y, + const std::vector& poly +) { + size_t n = poly.size(); + if (n < 3) return false; // no area → never inside + + bool inside = false; + for (size_t i = 0, j = n - 1; i < n; j = i++) { + double xi = poly[i].x, yi = poly[i].y; + double xj = poly[j].x, yj = poly[j].y; + + // does edge (j→i) straddle the horizontal ray at y? + bool cond = ((yi > y) != (yj > y)); + if (cond) { + std::cout<<"Checking point ("< 200) - { + // if (Verbosity() > 200) + // { std::cout << " populate phi bins for " << " layernum " << layernum << " phi " << phi + << " rad_gem " << rad_gem << " sigmaT " << sigmaT //<< " zigzag_pads " << zigzag_pads << std::endl; - } + //} std::vector pad_phibin; std::vector pad_phibin_share; populate_zigzag_phibins(side, layernum, phi, sigmaT, pad_phibin, pad_phibin_share); + + std::vector pad_phibin_serf; + std::vector pad_phibin_share_serf; + + SERF_zigzag_phibins(side, layernum, phi, rad_gem, sigmaT, pad_phibin_serf, pad_phibin_share_serf); /* if (pad_phibin.size() == 0) { */ /* pass_data.neff_electrons = 0; */ /* } else { */ @@ -897,39 +970,167 @@ void PHG4TpcPadPlaneReadout::build_serf_zigzag_phibins(const unsigned int side, } */ -/* -void getZigzagPadFractions( - double x0_mm, - double y0_mm, - const unsigned int layer, - double sigma_r, - std::map pad_fractions - ) { - double r0 = std::hypot(x0_mm, y0_mm); - int module = (layer - 7)/16; - if(module < 0) return {{}, module}; +double PHG4TpcPadPlaneReadout::integratedDensityOfCircleAndPad( + double hitX, + double hitY, + double sigma, + const std::vector& pad, + double gridStep +) { + // sanity checks + if (sigma <= 0.0 || pad.empty()) return 0.0; + + // constants + const int nSigma = 1; + double R = nSigma * sigma; + double gaussConst = 1.0 / (2.0 * M_PI * sigma * sigma); + double expDenominator = 2.0 * sigma * sigma; + + // pick a default grid step if none provided + if (gridStep <= 0.0) { + gridStep = sigma / 50.0; + } + + // 1) Compute pad bounding box + 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); + miny = std::min(miny, pad[i].y); + maxy = std::max(maxy, pad[i].y); + } + + // 2) Intersect that box with the 3σ circle’s box + double x0 = std::max(minx, hitX - R); + double x1 = std::min(maxx, hitX + R); + double y0 = std::max(miny, hitY - R); + double y1 = std::min(maxy, hitY + R); + + if (x1 <= x0 || y1 <= y0) { + // no overlap + return 0.0; + } + + // 3) Determine grid dimensions + int nx = static_cast(std::ceil((x1 - x0) / gridStep)); + int ny = static_cast(std::ceil((y1 - y0) / gridStep)); + double total = 0.0; + + // 4) Loop over cell centers + for (int ix = 0; ix < nx; ++ix) { + double x = x0 + (ix + 0.5) * gridStep; + for (int iy = 0; iy < ny; ++iy) { + double y = y0 + (iy + 0.5) * gridStep; + // skip outside the circle + double dx = x - hitX, dy = y - hitY; + if (dx*dx + dy*dy > R*R) continue; + // skip outside the pad polygon + if (!pointInPolygon(x, y, pad)) continue; + // accumulate Gaussian density + total += gaussConst * std::exp(-(dx*dx + dy*dy) / expDenominator); + } + } + + // 5) multiply by cell area for the approximate integral + return total * (gridStep * gridStep); +} + +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) +{ + const double radius = LayerGeom->get_radius(); + const double phistepsize = LayerGeom->get_phistep(); + const auto phibins = LayerGeom->get_phibins(); + + int sector1 = 0; + double x = rad_gem * cos(phi); + double y = rad_gem * sin(phi); + double xNew = 0.0; + double yNew = 0.0; + double phiNew = std::atan2(yNew,xNew); + double rad_gem_new = get_r(xNew, yNew); + rotatePointToSector( x, y, xNew, yNew, sector1); + std::cout<<"PHG4TpcPadPlaneReadout::populate_zigzag_phibins: x = "< 1e-8) pad_fractions[pad_id] = fraction; + + int phibin_low = LayerGeom->get_phibin(philim_high, side); + int phibin_high = LayerGeom->get_phibin(philim_low, side); + int npads = phibin_high - phibin_low; +/* + const double radlim_low_calc = rad_gem - (_nsigmas * cloud_sig_rp / radius) - phistepsize; + const double radlim_high_calc = rad_gem + (_nsigmas * cloud_sig_rp / radius) + phistepsize; + + // Find the pad range that covers this phi range + const double radlim_low = check_phi(side, phi, radlim_low_calc); + const double radlim_high = check_phi(side, phi, radlim_high_calc); + if (radlim_low < (LayerGeom->get_radius() - LayerGeom->get_thickness() / 2.0)) + { + std::cout << " radlim_low " << radlim_low << " is in the previous layer " <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; + + 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; } - //double total = 0; - //for(auto& kv : pad_fractions) total += kv.second; - //if(total>0) for(auto& kv : pad_fractions) kv.second /= total; - -}*/ + auto padinfo = Pads[layernum][pad_now - ntpc_phibins_sector[tpc_module]*sector]; + double charge = integratedDensityOfCircleAndPad( xNew, yNew, cloud_sig_rp / rad_gem_new, padinfo.vertices); + phibin_pad.push_back(pad_now); + phibin_pad_share.push_back(charge); + } + + 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) @@ -954,14 +1155,14 @@ void PHG4TpcPadPlaneReadout::populate_zigzag_phibins(const unsigned int side, co int phibin_high = LayerGeom->get_phibin(philim_low, side); int npads = phibin_high - phibin_low; - if (Verbosity() > 1000) - { - if (layernum == print_layer) - { + //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; - } - } + // } + // } if (npads < 0 || npads > 9) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h index 54aadd304a..26d2b4231c 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h @@ -55,7 +55,7 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane // 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 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 build_serf_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); double check_phi(const unsigned int side, const double phi, const double radius); @@ -123,17 +123,41 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane struct Point { double x, y; }; struct PadInfo { - std::string name; // pad name - int id; // pad id - double cx, cy; // centroid coords - std::vector vertices; // pad polygon + std::string name; // pad name + int pad_number; // pad number (number in module) + int pad_bin; // pad phi bin (number according to get_phi_bin) + double cx, cy; // centroid coords + double rad, phi; // pad radius and phi + std::vector vertices; // pad polygon + bool isedge = false; // whether to keep this pad signal + void clear() { + name.clear(); + pad_number = -1; + pad_bin = -1; + cx = cy = rad = phi = 0.0; + vertices.clear(); + isedge = false; + } }; -std::array,3> Pads; - +std::array,3*16+7> Pads; + double integratedDensityOfCircleAndPad(double hitX,double hitY, double sigma, const std::vector& pad,double gridStep = 0.0); // hard‑coded list of input .brd files static const std::vector brdMaps_; void loadPadPlanes(); +int ntpc_phibins_sector[3] = { 94, 128, 192 }; +bool pointInPolygon(double x, double y, std::vector poly); +int findPadForPoint( double x, double y, int tpc_module); + const std::array Thickness = + {{ + 0.56598621677629212, + 1.0206889851687158, + 1.0970475085472556, + 0.5630547309825637, + 0.56891770257002054, + }}; +double min_radii_module[3]={314.9836110818037, 416.59202613529567, 589.1096495597712}; +double max_radii_module[3]={399.85222874031024, 569.695373910603, 753.6667758418596}; }; From c50d8537d08f0e462d1df5dd46677445ab707c46 Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Fri, 8 Aug 2025 09:41:58 -0400 Subject: [PATCH 04/13] update pad shift --- .../g4tpc/PHG4TpcPadPlaneReadout.cc | 92 ++++++++++++++----- .../g4tpc/PHG4TpcPadPlaneReadout.h | 1 + 2 files changed, 69 insertions(+), 24 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index 8161e5f191..a6e1583535 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -205,11 +205,13 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) std::cout<<"Pad in module 0, pad 0: " << Pads[0][check_pad].name << " at (" << Pads[0][check_pad].cx << ", " << Pads[0][check_pad].cy << ")\n"; - for(size_t i = 0; i < Pads[0][check_pad].vertices.size(); i++) + for (int j=0; j<7+16*3;j++){ + for(size_t i = 0; i < Pads[j].size(); i++) { - std::cout<<"Vertex "<") == 0) { inSignal = false; keepSignal = false; - std::cout<<"Module "<& poly @@ -489,7 +491,7 @@ bool PHG4TpcPadPlaneReadout::pointInPolygon( } return inside; } -*/ + /* //------------------------------------------------------------------------------ // 2) findPadForPoint: returns padNumber or −1 if none @@ -728,6 +730,17 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( // Normalize the shares so they add up to 1 double norm1 = 0.0; + for (unsigned int ipad = 0; ipad < pad_phibin_serf.size(); ++ipad) + { + double pad_share = pad_phibin_share_serf[ipad]; + norm1 += pad_share; + } + for (unsigned int iphi = 0; iphi < pad_phibin_serf.size(); ++iphi) + { + pad_phibin_share_serf[iphi] /= norm1; + } + +norm1 = 0.0; for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad) { double pad_share = pad_phibin_share[ipad]; @@ -738,6 +751,12 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( pad_phibin_share[iphi] /= norm1; } + std::cout<<"--------------------------"< 100 && layernum == print_layer) @@ -981,7 +1000,7 @@ double PHG4TpcPadPlaneReadout::integratedDensityOfCircleAndPad( if (sigma <= 0.0 || pad.empty()) return 0.0; // constants - const int nSigma = 1; + const int nSigma = _nsigmas; double R = nSigma * sigma; double gaussConst = 1.0 / (2.0 * M_PI * sigma * sigma); double expDenominator = 2.0 * sigma * sigma; @@ -992,20 +1011,22 @@ double PHG4TpcPadPlaneReadout::integratedDensityOfCircleAndPad( } // 1) Compute pad bounding box - double minx = pad[0].x, maxx = pad[0].x; - double miny = pad[0].y, maxy = pad[0].y; + double minx = pad[0].x, maxx = pad.back().x; + double miny = pad[0].y, maxy = pad.back().y; for (size_t i = 1; i < pad.size(); ++i) { minx = std::min(minx, pad[i].x); maxx = std::max(maxx, pad[i].x); miny = std::min(miny, pad[i].y); maxy = std::max(maxy, pad[i].y); + // std::cout<<"Pad point "<get_phicenter(pad_now, side) <<" pad phi center "<get_phicenter(pad_keep[ipad], side) << std::endl; } return; diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h index 26d2b4231c..809b968754 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h @@ -141,6 +141,7 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane }; 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); // hard‑coded list of input .brd files static const std::vector brdMaps_; From 8ea43b56a22d57c0b294a38d160134fe51361b3f Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Fri, 15 Aug 2025 09:54:06 -0400 Subject: [PATCH 05/13] simulation/g4simulation/g4tpc/PHG4TpcDetector.cc --- .../g4simulation/g4tpc/PHG4TpcElectronDrift.cc | 10 +++++++++- .../g4simulation/g4tpc/PHG4TpcElectronDrift.h | 4 ++++ .../g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc | 13 +++++++++---- .../g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h | 4 +++- 4 files changed, 25 insertions(+), 6 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc index 617664100c..1bfaef4f88 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc @@ -355,6 +355,8 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) se->registerHisto(nthit); se->registerHisto(ntpad); } + h_adc_ref = new TH2F("h_adc_ref", "ADC Map. Coresoftware; Pad;Number of pads in cluster;ADC", 7, 0, 7.0, 7, 0.0, 7.0); + h_adc_serf= new TH2F("h_adc_serf", "ADC Map. SERF; Pad;Number of pads in cluster;ADC", 7, 0, 7.0, 7, 0.0, 7.0); padplane->InitRun(topNode); @@ -700,7 +702,7 @@ 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); + side, hiter, ntpad, nthit, h_adc_ref,h_adc_serf); } // end loop over electrons for this g4hit if (do_ElectronDriftQAHistos) @@ -926,6 +928,12 @@ int PHG4TpcElectronDrift::End(PHCompositeNode * /*topNode*/) ratioElectronsRR->Write(); EDrift_outf->Close(); } + + m_outf_adc.reset(new TFile("m_outf_adc.root", "recreate")); + m_outf_adc->cd(); + h_adc_ref->Write(); + h_adc_serf->Write(); + m_outf_adc->Close(); return Fun4AllReturnCodes::EVENT_OK; } diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h index 8bf3a3f7d4..6beaa85bb2 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h @@ -151,6 +151,10 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface std::unique_ptr m_outf; std::unique_ptr EDrift_outf; + std::unique_ptr m_outf_adc; + TH2 *h_adc_ref= nullptr; + TH2 *h_adc_serf= nullptr; + std::string detector; std::string hitnodename; std::string seggeonodename; diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index a6e1583535..9641d2a99b 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -212,6 +212,8 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) } } + + return Fun4AllReturnCodes::EVENT_OK; } @@ -515,7 +517,7 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( 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*/) + PHG4HitContainer::ConstIterator hiter, TNtuple * /*ntpad*/, TNtuple * /*nthit*/, TH2* h_adc_ref , TH2* h_adc_serf ) { // One electron per call of this method // The x_gem and y_gem values have already been randomized within the transverse drift diffusion width @@ -754,6 +756,8 @@ norm1 = 0.0; std::cout<<"--------------------------"<Fill(iphi, pad_phibin.size(), pad_phibin_share[iphi]); + h_adc_serf->Fill(iphi, pad_phibin_serf.size(), pad_phibin_share_serf[iphi]); std::cout<<" phibin ref "< Date: Fri, 15 Aug 2025 09:56:48 -0400 Subject: [PATCH 06/13] Added histos --- simulation/g4simulation/g4tpc/PHG4TpcDetector.cc | 5 +++-- simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc b/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc index 93819148e2..04573e6e53 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcDetector.cc @@ -537,8 +537,9 @@ void PHG4TpcDetector::add_geometry_node() } m_cdb = CDBInterface::instance(); - std::string calibdir = m_cdb->getUrl("TPC_FEE_CHANNEL_MAP"); - + //std::string calibdir = m_cdb->getUrl("TPC_FEE_CHANNEL_MAP"); + std::string calibdir = "/sphenix/user/mitrankova/PadPlane_Readout/map/TPCPadPlaneCDBTTree.root"; + if (! calibdir.empty()) { m_cdbttree = new CDBTTree(calibdir); diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h index 4d2a9371e6..f6de5e99df 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h @@ -16,6 +16,7 @@ class TrkrHitTruthAssoc; class PHCompositeNode; class PHG4TpcCylinderGeomContainer; class TNtuple; +class TH2; class PHG4TpcPadPlane : public SubsysReco, public PHParameterInterface { @@ -30,7 +31,7 @@ class PHG4TpcPadPlane : public SubsysReco, public PHParameterInterface int InitRun(PHCompositeNode *topNode) override; virtual void UpdateInternalParameters() { return; } // virtual void MapToPadPlane(PHG4CellContainer * /*g4cells*/, const double /*x_gem*/, const double /*y_gem*/, const double /*t_gem*/, const unsigned int /*side*/, PHG4HitContainer::ConstIterator /*hiter*/, TNtuple * /*ntpad*/, TNtuple * /*nthit*/) {} - virtual void MapToPadPlane(TpcClusterBuilder& /*builder*/, 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*/)=0;// { return {}; } + virtual void MapToPadPlane(TpcClusterBuilder& /*builder*/, 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 )=0;// { return {}; } void Detector(const std::string &name) { detector = name; } protected: From 0edba28707e866360c700a14a1965d8ab642576b Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Sat, 23 Aug 2025 23:07:59 -0400 Subject: [PATCH 07/13] Works, but r is not taken into account --- .../g4tpc/PHG4TpcElectronDrift.cc | 16 ++- .../g4simulation/g4tpc/PHG4TpcElectronDrift.h | 6 +- .../g4simulation/g4tpc/PHG4TpcPadPlane.h | 5 +- .../g4tpc/PHG4TpcPadPlaneReadout.cc | 134 ++++++++++-------- .../g4tpc/PHG4TpcPadPlaneReadout.h | 3 +- 5 files changed, 97 insertions(+), 67 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc index 1bfaef4f88..b66e43952c 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.cc @@ -355,8 +355,8 @@ int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode) se->registerHisto(nthit); se->registerHisto(ntpad); } - h_adc_ref = new TH2F("h_adc_ref", "ADC Map. Coresoftware; Pad;Number of pads in cluster;ADC", 7, 0, 7.0, 7, 0.0, 7.0); - h_adc_serf= new TH2F("h_adc_serf", "ADC Map. SERF; Pad;Number of pads in cluster;ADC", 7, 0, 7.0, 7, 0.0, 7.0); + //h_adc_ref = new TH2F("h_adc_ref", "ADC Map. Coresoftware; Pad;Number of pads in cluster;ADC", 7, 0, 7.0, 7, 0.0, 7.0); + //h_adc_serf= new TH2F("h_adc_serf", "ADC Map. SERF; Pad;Number of pads in cluster;ADC", 7, 0, 7.0, 7, 0.0, 7.0); padplane->InitRun(topNode); @@ -701,8 +701,13 @@ int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode) nt->Fill(ihit, t_start, t_final, t_sigma, rad_final, z_start, z_final); } padplane->MapToPadPlane(truth_clusterer, single_hitsetcontainer.get(), + temp_hitsetcontainer.get(), hittruthassoc, x_final, y_final, t_final, + side, hiter, ntpad, nthit); + /*padplane->MapToPadPlane(truth_clusterer, single_hitsetcontainer.get(), temp_hitsetcontainer.get(), hittruthassoc, x_final, y_final, t_final, side, hiter, ntpad, nthit, h_adc_ref,h_adc_serf); + */ + // std::cout<<" Electron Drift h_adc_ref->GetEntries() = "<GetEntries()<<" h_adc_serf->GetEntries() = "<GetEntries()<Close(); } - m_outf_adc.reset(new TFile("m_outf_adc.root", "recreate")); + /*m_outf_adc.reset(new TFile("m_outf_adc.root", "recreate")); + if (m_outf_adc) + std::cout<<" file m_outf_adc.root created "<cd(); h_adc_ref->Write(); h_adc_serf->Write(); m_outf_adc->Close(); + */ return Fun4AllReturnCodes::EVENT_OK; } diff --git a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h index 6beaa85bb2..004c32001b 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcElectronDrift.h @@ -151,9 +151,9 @@ class PHG4TpcElectronDrift : public SubsysReco, public PHParameterInterface std::unique_ptr m_outf; std::unique_ptr EDrift_outf; - std::unique_ptr m_outf_adc; - TH2 *h_adc_ref= nullptr; - TH2 *h_adc_serf= nullptr; + //std::unique_ptr m_outf_adc; + //TH2 *h_adc_ref= nullptr; + //TH2 *h_adc_serf= nullptr; std::string detector; std::string hitnodename; diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h index f6de5e99df..196b54797e 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlane.h @@ -16,7 +16,7 @@ class TrkrHitTruthAssoc; class PHCompositeNode; class PHG4TpcCylinderGeomContainer; class TNtuple; -class TH2; +//class TH2; class PHG4TpcPadPlane : public SubsysReco, public PHParameterInterface { @@ -31,7 +31,8 @@ class PHG4TpcPadPlane : public SubsysReco, public PHParameterInterface int InitRun(PHCompositeNode *topNode) override; virtual void UpdateInternalParameters() { return; } // virtual void MapToPadPlane(PHG4CellContainer * /*g4cells*/, const double /*x_gem*/, const double /*y_gem*/, const double /*t_gem*/, const unsigned int /*side*/, PHG4HitContainer::ConstIterator /*hiter*/, TNtuple * /*ntpad*/, TNtuple * /*nthit*/) {} - virtual void MapToPadPlane(TpcClusterBuilder& /*builder*/, 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 )=0;// { return {}; } + virtual void MapToPadPlane(TpcClusterBuilder& /*builder*/, 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*/)=0;// { return {}; } + //virtual void MapToPadPlane(TpcClusterBuilder& /*builder*/, 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 )=0;// { return {}; } void Detector(const std::string &name) { detector = name; } protected: diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index 9641d2a99b..ce78ba0b52 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -196,16 +196,9 @@ int PHG4TpcPadPlaneReadout::InitRun(PHCompositeNode *topNode) // 3) (optional) print a summary - for (size_t i = 0; i < Pads.size(); ++i) - { - std::cout << "[PHG4TpcPadPlaneReadout] Module " << i - << " has " << Pads[i].size() << " pads\n"; - } - int check_pad = 81; - std::cout<<"Pad in module 0, pad 0: " - << Pads[0][check_pad].name << " at (" - << Pads[0][check_pad].cx << ", " << Pads[0][check_pad].cy << ")\n"; - for (int j=0; j<7+16*3;j++){ + + + /* 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 = "<") == 0) { inSignal = false; keepSignal = false; - //std::cout<<"Module "< pad_phibin; std::vector pad_phibin_share; - populate_zigzag_phibins(side, layernum, phi, sigmaT, pad_phibin, pad_phibin_share); + //std::vector pad_phibin_ref; + //std::vector pad_phibin_share_ref; + //std::vector pad_phibin_serf; + //std::vector pad_phibin_share_serf; - std::vector pad_phibin_serf; - std::vector pad_phibin_share_serf; + // populate_zigzag_phibins(side, layernum, phi, sigmaT, pad_phibin_ref, pad_phibin_share_ref); - SERF_zigzag_phibins(side, layernum, phi, rad_gem, sigmaT, pad_phibin_serf, pad_phibin_share_serf); + + populate_zigzag_phibins(side, layernum, phi, sigmaT, pad_phibin, pad_phibin_share); + + //SERF_zigzag_phibins(side, layernum, phi, rad_gem, sigmaT, pad_phibin, pad_phibin_share); /* if (pad_phibin.size() == 0) { */ /* pass_data.neff_electrons = 0; */ /* } else { */ @@ -732,35 +732,35 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( // Normalize the shares so they add up to 1 double norm1 = 0.0; - for (unsigned int ipad = 0; ipad < pad_phibin_serf.size(); ++ipad) + for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad) { - double pad_share = pad_phibin_share_serf[ipad]; + double pad_share = pad_phibin_share[ipad]; norm1 += pad_share; } - for (unsigned int iphi = 0; iphi < pad_phibin_serf.size(); ++iphi) + for (unsigned int iphi = 0; iphi < pad_phibin.size(); ++iphi) { - pad_phibin_share_serf[iphi] /= norm1; + pad_phibin_share[iphi] /= norm1; } - +/* norm1 = 0.0; - for (unsigned int ipad = 0; ipad < pad_phibin.size(); ++ipad) + 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<<"--------------------------"<Fill(iphi, pad_phibin.size(), pad_phibin_share[iphi]); - h_adc_serf->Fill(iphi, pad_phibin_serf.size(), pad_phibin_share_serf[iphi]); - std::cout<<" phibin ref "< 100 && layernum == print_layer) @@ -1079,11 +1079,11 @@ void PHG4TpcPadPlaneReadout::SERF_zigzag_phibins(const unsigned int side, const double yNew = 0.0; rotatePointToSector( x, y, xNew, yNew, sector1); - 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; } */ - 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; + // 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; for (int ipad = 0; ipad <= npads; ipad++) { @@ -1150,10 +1149,17 @@ std::cout<<"!!!!!! phi = "<= ntpc_phibins_sector[tpc_module]){ + look_pad-=ntpc_phibins_sector[tpc_module]; + } + + look_pad = ntpc_phibins_sector[tpc_module] - look_pad -1; + + // std::cout<<"pad now = "<get_phicenter(pad_now, side) <<" pad phi center "<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) - //{ + if (Verbosity() > 1000) + { + if (layernum == print_layer) + { std::cout << " REF zigzags: phi " << phi << " philim_low " << philim_low << " phibin_low " << phibin_low - << " philim_high " << philim_high << " phibin_high " << phibin_high << " npads " << npads << std::endl; - // } - // } + << " philim_high " << philim_high << " phibin_high " << phibin_high << " npads " << npads <<" sector "< 9) { diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h index a8bc2adbcc..bf4860c098 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.h @@ -45,7 +45,8 @@ class PHG4TpcPadPlaneReadout : public PHG4TpcPadPlane // otherwise warning of inconsistent overload since only one MapToPadPlane methow is overridden using PHG4TpcPadPlane::MapToPadPlane; - 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 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; From 9221806017953a2f5bab0527eb7ca7e231b80735 Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Sat, 23 Aug 2025 23:08:24 -0400 Subject: [PATCH 08/13] correct side typo --- simulation/g4simulation/g4eval/SvtxTruthEval.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simulation/g4simulation/g4eval/SvtxTruthEval.cc b/simulation/g4simulation/g4eval/SvtxTruthEval.cc index f36135399b..e24f672852 100644 --- a/simulation/g4simulation/g4eval/SvtxTruthEval.cc +++ b/simulation/g4simulation/g4eval/SvtxTruthEval.cc @@ -737,7 +737,7 @@ void SvtxTruthEval::G4ClusterSize(TrkrDefs::cluskey ckey, unsigned int layer, co double avge_z = (outer_z + inner_z) / 2.0; unsigned int side = 0; - if (avge_z < 0) + if (avge_z > 0) { side = 1; } From acb10e39b4d186a164bd13bedcefd4e20587e092 Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Sat, 23 Aug 2025 23:09:12 -0400 Subject: [PATCH 09/13] Added check printouts --- offline/packages/tpc/TpcClusterizer.cc | 35 +++++++++++- .../tpc/TpcCombinedRawDataUnpacker.cc | 53 +++++++++++++++++-- .../packages/tpc/TpcCombinedRawDataUnpacker.h | 30 +++++++++-- offline/packages/trackreco/PHCASeeding.cc | 1 + 4 files changed, 112 insertions(+), 7 deletions(-) diff --git a/offline/packages/tpc/TpcClusterizer.cc b/offline/packages/tpc/TpcClusterizer.cc index 7b9110448e..b3a6e6785b 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -683,6 +683,39 @@ namespace // std::cout << "done transform" << std::endl; // we need the cluster key and all associated hit keys (note: the cluster key includes the hitset key) +/* +pthread_mutex_lock(&mythreadlock); + +std::cout << "Cluster side = " << my_data.side + << " layer = " << my_data.layer + << " sector = " << my_data.sector + << " tpcHitSetKey = " << tpcHitSetKey + << " subsurfkey = " << subsurfkey + << " clusx = " << clusx + << " clusy = " << clusy + << " clusz = " << clusz + << " local X = " << local(0) + << " local Y = " << clust + << " center phi = " << clusphi + << ", made from hits: " << std::endl; + +for (auto const& hit : ihit_list) +{ + double hitphi = my_data.layergeom->get_phi(hit.iphi + my_data.phioffset); + std::cout << hitphi << " "; +} +std::cout << std::endl; + +std::cout << "Diff: clusphi - hitphi = " << std::endl; +for (auto const& hit : ihit_list) +{ + double diff = clusphi - my_data.layergeom->get_phi(hit.iphi + my_data.phioffset); + std::cout << diff << " "; +} +std::cout << std::endl; + +pthread_mutex_unlock(&mythreadlock); +*/ TrkrCluster *clus_base = nullptr; bool b_made_cluster{false}; @@ -1654,7 +1687,7 @@ int TpcClusterizer::process_event(PHCompositeNode *topNode) // insert in map // std::cout << "X: " << cluster->getLocalX() << "Y: " << cluster->getLocalY() << std::endl; m_clusterlist->addClusterSpecifyKey(ckey, cluster); - +//std::cout<< "TpcClusterizer::process_event: added cluster : side = "<< data.side <<" sector = "<Branch("event", &t_event, "t_event/I"); + m_ntup_hits_tree->Branch("gtmbco", &t_gtmbco, "t_gtmbco/l"); // unsigned 64-bit + m_ntup_hits_tree->Branch("packid", &t_packid, "t_packid/i"); + m_ntup_hits_tree->Branch("sector", &t_sector, "t_sector/I"); + m_ntup_hits_tree->Branch("side", &t_side, "t_side/I"); + m_ntup_hits_tree->Branch("layer", &t_layer, "t_layer/I"); + m_ntup_hits_tree->Branch("fee", &t_fee, "t_fee/I"); + m_ntup_hits_tree->Branch("chan", &t_chan, "t_chan/I"); + m_ntup_hits_tree->Branch("sampadd", &t_sampadd, "t_sampadd/s"); // unsigned 16-bit + m_ntup_hits_tree->Branch("sampch", &t_sampch, "t_sampch/s"); + m_ntup_hits_tree->Branch("phi", &t_phi, "t_phi/F"); + m_ntup_hits_tree->Branch("phi_center", &t_phi_center,"t_phi_center/F"); + m_ntup_hits_tree->Branch("phibin", &t_phibin, "t_phibin/I"); + m_ntup_hits_tree->Branch("tbin", &t_tbin, "t_tbin/I"); + m_ntup_hits_tree->Branch("hitkey", &t_hitkey, "t_hitkey/i"); // unsigned 32-bit + m_ntup_hits_tree->Branch("adc", &t_adc, "t_adc/F"); + m_ntup_hits_tree->Branch("ped", &t_ped, "t_ped/F"); + m_ntup_hits_tree->Branch("width", &t_width, "t_width/F"); + + if (m_doChanHitsCut) { m_HitChanDis = new TH2F("HitChanDis", "HitChanDis", 451, -0.5, 450.5, 256, -0.5, 255.5); @@ -279,6 +301,7 @@ int TpcCombinedRawDataUnpacker::process_event(PHCompositeNode* topNode) double phi = ((side == 1 ? 1 : -1) * (m_cdbttree->GetDoubleValue(key, varname) - M_PI / 2.)) + ((sector % 12) * M_PI / 6); PHG4TpcCylinderGeom* layergeom = geom_container->GetLayerCellGeom(layer); unsigned int phibin = layergeom->get_phibin(phi, side); + //std::cout<<"TpcCombinedRawDataUnpacker:: side = "<::quiet_NaN(); + int t_sector= std::numeric_limits::quiet_NaN(); + int t_side= std::numeric_limits::quiet_NaN(); + int t_layer= std::numeric_limits::quiet_NaN(); + int t_fee= std::numeric_limits::quiet_NaN(); + int t_chan= std::numeric_limits::quiet_NaN(); + int t_phibin= std::numeric_limits::quiet_NaN(); + int t_tbin= std::numeric_limits::quiet_NaN(); + float t_phi=std::numeric_limits::quiet_NaN(); + float t_phi_center=std::numeric_limits::quiet_NaN(); + float t_adc=std::numeric_limits::quiet_NaN(); + float t_ped=std::numeric_limits::quiet_NaN(); + float t_width=std::numeric_limits::quiet_NaN(); + + uint32_t t_hitkey=std::numeric_limits::quiet_NaN(); + uint64_t t_gtmbco=std::numeric_limits::quiet_NaN(); + int32_t t_packid=std::numeric_limits::quiet_NaN(); + uint16_t t_sampadd=std::numeric_limits::quiet_NaN(); + uint16_t t_sampch=std::numeric_limits::quiet_NaN(); + + + bool m_writeTree{true}; bool m_do_baseline_corr{false}; int m_baseline_nsigma{2}; bool m_do_zs_emulation{false}; diff --git a/offline/packages/trackreco/PHCASeeding.cc b/offline/packages/trackreco/PHCASeeding.cc index 63030a86c7..09506da1eb 100644 --- a/offline/packages/trackreco/PHCASeeding.cc +++ b/offline/packages/trackreco/PHCASeeding.cc @@ -275,6 +275,7 @@ std::pair PHCASeeding::F // get global position, convert to Acts::Vector3 and store in map const Acts::Vector3 globalpos_d = getGlobalPosition(ckey, cluster); const Acts::Vector3 globalpos = {globalpos_d.x(), globalpos_d.y(), globalpos_d.z()}; + // std::cout<<"PHCASeeding::FillGlobalPositions: side: "<getZSize()<<" layer = "< Date: Sat, 23 Aug 2025 23:09:33 -0400 Subject: [PATCH 10/13] Chenxi output --- .../TrackingDiagnostics/TrackResiduals.cc | 68 ++++++++++++++++++- .../TrackingDiagnostics/TrackResiduals.h | 15 +++- 2 files changed, 78 insertions(+), 5 deletions(-) diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.cc b/offline/packages/TrackingDiagnostics/TrackResiduals.cc index ca036914aa..fed1387613 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.cc +++ b/offline/packages/TrackingDiagnostics/TrackResiduals.cc @@ -15,6 +15,8 @@ #include #include #include +#include +#include #include #include @@ -194,6 +196,8 @@ void TrackResiduals::clearClusterStateVectors() m_clusAdc.clear(); m_clusMaxAdc.clear(); m_cluslayer.clear(); + m_clussize.clear(); + m_clushitsetkey.clear(); m_statelx.clear(); m_statelz.clear(); @@ -219,6 +223,22 @@ int TrackResiduals::process_event(PHCompositeNode* topNode) auto *mvtxGeom = findNode::getClass(topNode, "CYLINDERGEOM_MVTX"); auto *inttGeom = findNode::getClass(topNode, "CYLINDERGEOM_INTT"); auto *mmGeom = findNode::getClass(topNode, "CYLINDERGEOM_MICROMEGAS_FULL"); + auto clusterhitassocmap = findNode::getClass(topNode, "TRKR_CLUSTERHITASSOC"); + auto clustercrossingassoc = findNode::getClass(topNode, "TRKR_CLUSTERCROSSINGASSOC"); + if (!clustercrossingassoc) + { + std::cerr << "ERROR: Can't find TRKR_CLUSTERCROSSINGASSOC node!" << std::endl; + return Fun4AllReturnCodes::ABORTEVENT; + } + if (!clusterhitassocmap) + { + std::cerr << "ERROR: Can't find TRKR_CLUSTERHITASSOC node!" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + std::cout << "TRKR_CLUSTERHITASSOC size: " << clusterhitassocmap->size() << std::endl; + + m_vx = m_vy = m_vz = std::numeric_limits::quiet_NaN(); if (!mmGeom) { @@ -310,7 +330,7 @@ int TrackResiduals::process_event(PHCompositeNode* topNode) if (m_doClusters) { - fillClusterTree(clustermap, geometry); + fillClusterTree(clusterhitassocmap, clustermap, clustercrossingassoc, geometry); } if (m_convertSeeds) @@ -675,13 +695,18 @@ void TrackResiduals::lineFitClusters(std::vector& keys, m_yzslope = std::get<0>(yzparams); } -void TrackResiduals::fillClusterTree(TrkrClusterContainer* clusters, - ActsGeometry* geometry) +void TrackResiduals::fillClusterTree(TrkrClusterHitAssoc *clusterhitassoc, TrkrClusterContainer *clusters, + TrkrClusterCrossingAssoc *clustercrossingassoc, ActsGeometry *geometry) { if (clusters->size() < m_min_cluster_size) { return; } + if (!clusterhitassoc) + { + std::cerr << "ERROR: Can't find TRKR_CLUSTERHITASSOC node!" << std::endl; + return; + } for (const auto& det : {TrkrDefs::TrkrId::mvtxId, TrkrDefs::TrkrId::inttId, TrkrDefs::TrkrId::tpcId, TrkrDefs::TrkrId::micromegasId}) { @@ -689,11 +714,14 @@ void TrackResiduals::fillClusterTree(TrkrClusterContainer* clusters, { m_scluslayer = TrkrDefs::getLayer(hitsetkey); auto range = clusters->getClusters(hitsetkey); + m_clustHitsetkey = std::numeric_limits::max(); + m_clustHitsetkey = hitsetkey; for (auto iter = range.first; iter != range.second; ++iter) { auto key = iter->first; auto *cluster = clusters->findCluster(key); Acts::Vector3 glob; + m_scluskey = key; // NOT IMPLEMENTED YET // if (TrkrDefs::getTrkrId(key) == TrkrDefs::tpcId) // { @@ -717,9 +745,17 @@ void TrackResiduals::fillClusterTree(TrkrClusterContainer* clusters, auto para_errors = m_clusErrPara.get_clusterv5_modified_error(cluster, m_sclusgr, key); m_phisize = cluster->getPhiSize(); m_zsize = cluster->getZSize(); + m_aclussize = cluster->getSize(); m_scluselx = std::sqrt(para_errors.first); m_scluselz = std::sqrt(para_errors.second); + m_clust_crossings.clear(); + auto crossingRange = clustercrossingassoc->getCrossings(key); + for (auto cxit = crossingRange.first; cxit != crossingRange.second; ++cxit) + { + short int crossing_number = cxit->second; + m_clust_crossings.push_back(crossing_number); + } //! Fill relevant geom info that is specific to subsystem switch (det) { @@ -778,7 +814,17 @@ void TrackResiduals::fillClusterTree(TrkrClusterContainer* clusters, default: break; } + m_clust_hitkeys.clear(); + auto hitrange = clusterhitassoc->getHits(key); + // size_t num_hits = std::distance(hitrange.first, hitrange.second); + // std::cout << "Cluster Key: " << key << " has " << num_hits << " associated hits." << std::endl; + for (auto hit_iter = hitrange.first; hit_iter != hitrange.second; ++hit_iter) + { + TrkrDefs::hitkey hkey = hit_iter->second; + m_clust_hitkeys.push_back(hkey); + // std::cout << " Hitkey: " << hkey << std::endl; + } m_clustree->Fill(); } } @@ -821,6 +867,7 @@ void TrackResiduals::fillHitTree(TrkrHitSetContainer* hitmap, PHG4CylinderGeomContainer* inttGeom, PHG4CylinderGeomContainer* mmGeom) { + m_hitkey = std::numeric_limits::max(); if (!tpcGeom or !mvtxGeom or !inttGeom or !mmGeom) { std::cout << PHWHERE << "missing hit map, can't continue with hit tree" @@ -967,6 +1014,7 @@ void TrackResiduals::fillHitTree(TrkrHitSetContainer* hitmap, } case TrkrDefs::TrkrId::tpcId: { + m_hitkey = hitkey; m_row = std::numeric_limits::quiet_NaN(); m_col = std::numeric_limits::quiet_NaN(); m_segtype = std::numeric_limits::quiet_NaN(); @@ -984,6 +1032,7 @@ void TrackResiduals::fillHitTree(TrkrHitSetContainer* hitmap, m_hitgx = glob.x(); m_hitgy = glob.y(); m_hitgz = glob.z(); + m_hit_phi = phi; break; } @@ -1209,6 +1258,8 @@ void TrackResiduals::fillClusterBranchesKF(TrkrDefs::cluskey ckey, SvtxTrack* tr m_cluslayer.push_back(TrkrDefs::getLayer(ckey)); m_clusphisize.push_back(cluster->getPhiSize()); m_cluszsize.push_back(cluster->getZSize()); + m_clussize.push_back(cluster->getPhiSize() * cluster->getZSize()); + m_clushitsetkey.push_back(TrkrDefs::getHitSetKeyFromClusKey(ckey)); auto misaligncenter = surf->center(geometry->geometry().getGeoContext()); auto misalignnorm = -1 * surf->normal(geometry->geometry().getGeoContext(), Acts::Vector3(1, 1, 1), Acts::Vector3(1, 1, 1)); @@ -1455,6 +1506,8 @@ void TrackResiduals::fillClusterBranchesSeeds(TrkrDefs::cluskey ckey, // SvtxTr m_cluslayer.push_back(TrkrDefs::getLayer(ckey)); m_clusphisize.push_back(cluster->getPhiSize()); m_cluszsize.push_back(cluster->getZSize()); + m_clussize.push_back(cluster->getPhiSize() * cluster->getZSize()); + m_clushitsetkey.push_back(TrkrDefs::getHitSetKeyFromClusKey(ckey)); if (Verbosity() > 1) { @@ -1690,9 +1743,11 @@ void TrackResiduals::createBranches() m_hittree->Branch("event", &m_event, "m_event/I"); m_hittree->Branch("gl1bco", &m_bco, "m_bco/l"); m_hittree->Branch("hitsetkey", &m_hitsetkey, "m_hitsetkey/i"); + m_hittree->Branch("hitkey", &m_hitkey, "m_hitkey/i"); m_hittree->Branch("gx", &m_hitgx, "m_hitgx/F"); m_hittree->Branch("gy", &m_hitgy, "m_hitgy/F"); m_hittree->Branch("gz", &m_hitgz, "m_hitgz/F"); + m_hittree->Branch("phi", &m_hit_phi, "m_hit_phi/F"); m_hittree->Branch("layer", &m_hitlayer, "m_hitlayer/I"); m_hittree->Branch("sector", &m_sector, "m_sector/I"); m_hittree->Branch("side", &m_side, "m_side/I"); @@ -1714,6 +1769,9 @@ void TrackResiduals::createBranches() m_hittree->Branch("mbdcharge",&m_totalmbd, "m_totalmbd/F"); m_clustree = new TTree("clustertree", "A tree with all clusters"); + m_clustree->Branch("cluskey", &m_scluskey); + m_clustree->Branch("hitsetkey", &m_clustHitsetkey, "m_clustHitsetkey/i"); + m_clustree->Branch("clus_hitkeys", &m_clust_hitkeys); m_clustree->Branch("run", &m_runnumber, "m_runnumber/I"); m_clustree->Branch("segment", &m_segment, "m_segment/I"); m_clustree->Branch("event", &m_event, "m_event/I"); @@ -1722,9 +1780,12 @@ void TrackResiduals::createBranches() m_clustree->Branch("lz", &m_scluslz, "m_scluslz/F"); m_clustree->Branch("gx", &m_sclusgx, "m_sclusgx/F"); m_clustree->Branch("gy", &m_sclusgy, "m_sclusgy/F"); + m_clustree->Branch("r", &m_sclusgr, "m_sclusgr/F"); m_clustree->Branch("gz", &m_sclusgz, "m_sclusgz/F"); m_clustree->Branch("phi", &m_sclusphi, "m_sclusphi/F"); m_clustree->Branch("eta", &m_scluseta, "m_scluseta/F"); + m_clustree->Branch("clussize", &m_aclussize, "m_aclussize/I"); + m_clustree->Branch("layer", &m_scluslayer, "m_scluslayer/I"); m_clustree->Branch("adc", &m_adc, "m_adc/F"); m_clustree->Branch("phisize", &m_phisize, "m_phisize/I"); m_clustree->Branch("zsize", &m_zsize, "m_zsize/I"); @@ -1741,6 +1802,7 @@ void TrackResiduals::createBranches() m_clustree->Branch("timebucket", &m_timebucket, "m_timebucket/I"); m_clustree->Branch("segtype", &m_segtype, "m_segtype/I"); m_clustree->Branch("tile", &m_tileid, "m_tileid/I"); + m_clustree->Branch("clust_crossings", &m_clust_crossings); m_tree = new TTree("residualtree", "A tree with track, cluster, and state info"); m_tree->Branch("run", &m_runnumber, "m_runnumber/I"); diff --git a/offline/packages/TrackingDiagnostics/TrackResiduals.h b/offline/packages/TrackingDiagnostics/TrackResiduals.h index 983a3f5b4b..10ac6c6fa8 100644 --- a/offline/packages/TrackingDiagnostics/TrackResiduals.h +++ b/offline/packages/TrackingDiagnostics/TrackResiduals.h @@ -30,6 +30,8 @@ class TrkrHitSetContainer; class PHG4TpcCylinderGeomContainer; class PHG4CylinderGeomContainer; class TpcDistortionCorrectionContainer; +class TrkrClusterHitAssoc; +class TrkrClusterCrossingAssoc; class TrackResiduals : public SubsysReco { public: @@ -72,7 +74,7 @@ class TrackResiduals : public SubsysReco void createBranches(); static float convertTimeToZ(ActsGeometry *geometry, TrkrDefs::cluskey cluster_key, TrkrCluster *cluster); void fillEventTree(PHCompositeNode *topNode); - void fillClusterTree(TrkrClusterContainer *clusters, ActsGeometry *geometry); + void fillClusterTree(TrkrClusterHitAssoc *clusterhitassoc, TrkrClusterContainer *clusters, TrkrClusterCrossingAssoc *clustercrossingassoc, ActsGeometry *geometry); void fillHitTree(TrkrHitSetContainer *hitmap, ActsGeometry *geometry, PHG4TpcCylinderGeomContainer *tpcGeom, PHG4CylinderGeomContainer *mvtxGeom, PHG4CylinderGeomContainer *inttGeom, PHG4CylinderGeomContainer *mmGeom); @@ -221,9 +223,11 @@ class TrackResiduals : public SubsysReco //! hit tree info uint32_t m_hitsetkey = std::numeric_limits::quiet_NaN(); + TrkrDefs::hitkey m_hitkey = std::numeric_limits::max(); float m_hitgx = std::numeric_limits::quiet_NaN(); float m_hitgy = std::numeric_limits::quiet_NaN(); float m_hitgz = std::numeric_limits::quiet_NaN(); + float m_hit_phi = std::numeric_limits::quiet_NaN(); int m_hitlayer = std::numeric_limits::quiet_NaN(); int m_sector = std::numeric_limits::quiet_NaN(); int m_hitpad = std::numeric_limits::quiet_NaN(); @@ -239,6 +243,9 @@ class TrackResiduals : public SubsysReco int m_nvertices = std::numeric_limits::quiet_NaN(); //! cluster tree info + TrkrDefs::cluskey m_scluskey; + std::vector m_clust_hitkeys; + uint32_t m_clustHitsetkey = std::numeric_limits::max(); float m_sclusgr = std::numeric_limits::quiet_NaN(); float m_sclusphi = std::numeric_limits::quiet_NaN(); float m_scluseta = std::numeric_limits::quiet_NaN(); @@ -246,6 +253,7 @@ class TrackResiduals : public SubsysReco float m_clusmaxadc = std::numeric_limits::quiet_NaN(); int m_phisize = std::numeric_limits::quiet_NaN(); int m_zsize = std::numeric_limits::quiet_NaN(); + int m_aclussize = std::numeric_limits::quiet_NaN(); float m_scluslx = std::numeric_limits::quiet_NaN(); float m_scluslz = std::numeric_limits::quiet_NaN(); float m_sclusgx = std::numeric_limits::quiet_NaN(); @@ -282,11 +290,13 @@ class TrackResiduals : public SubsysReco std::vector m_clsector; std::vector m_clside; std::vector m_cluslayer; + std::vector m_clussize; std::vector m_clusphisize; std::vector m_cluszsize; std::vector m_clusedge; std::vector m_clusoverlap; std::vector m_cluskeys; + std::vector m_clushitsetkey; std::vector m_idealsurfcenterx; std::vector m_idealsurfcentery; std::vector m_idealsurfcenterz; @@ -308,7 +318,8 @@ class TrackResiduals : public SubsysReco std::vector m_missurfalpha; std::vector m_missurfbeta; std::vector m_missurfgamma; - + std::vector m_clust_crossings; + //! states on track information std::vector m_statelx; std::vector m_statelz; From 63866581f8909d0a4472f1bd101453d1036d0955 Mon Sep 17 00:00:00 2001 From: Mariia Mitrankova Date: Tue, 26 Aug 2025 21:05:07 -0400 Subject: [PATCH 11/13] more printout --- offline/packages/tpc/TpcClusterizer.cc | 43 ++++++++++++++++++++++---- 1 file changed, 37 insertions(+), 6 deletions(-) diff --git a/offline/packages/tpc/TpcClusterizer.cc b/offline/packages/tpc/TpcClusterizer.cc index b3a6e6785b..b488617013 100644 --- a/offline/packages/tpc/TpcClusterizer.cc +++ b/offline/packages/tpc/TpcClusterizer.cc @@ -683,9 +683,9 @@ namespace // std::cout << "done transform" << std::endl; // we need the cluster key and all associated hit keys (note: the cluster key includes the hitset key) -/* -pthread_mutex_lock(&mythreadlock); +pthread_mutex_lock(&mythreadlock); +std::cout << "==============================================" << std::endl; std::cout << "Cluster side = " << my_data.side << " layer = " << my_data.layer << " sector = " << my_data.sector @@ -696,9 +696,19 @@ std::cout << "Cluster side = " << my_data.side << " clusz = " << clusz << " local X = " << local(0) << " local Y = " << clust + << " cluster ADC = " << adc_sum + << " center iphibin = " << clusiphi << " center phi = " << clusphi - << ", made from hits: " << std::endl; - + << " " << std::endl; +std::cout << "-----------------------------------------------" << std::endl; +std::cout << " Sector min phi = " << my_data.layergeom->get_sector_min_phi()[my_data.side][my_data.sector] + <<" Sector max phi = " << my_data.layergeom->get_sector_max_phi()[my_data.side][my_data.sector] + <<" Sector bin width = " << my_data.layergeom->get_phistep() + <<" nphibins = "<get_phibins() + << std::endl; +std::cout << "-----------------------------------------------" << std::endl; + +std::cout << "Hit pad bins = " << std::endl; for (auto const& hit : ihit_list) { double hitphi = my_data.layergeom->get_phi(hit.iphi + my_data.phioffset); @@ -706,6 +716,27 @@ for (auto const& hit : ihit_list) } std::cout << std::endl; +std::cout << "Hit pad bins = " << std::endl; +for (auto const& hit : ihit_list) +{ + std::cout << hit.iphi + my_data.phioffset << " "; +} +std::cout << std::endl; + +std::cout << "Hitkeys = " << std::endl; +for (size_t i=0; i Date: Tue, 26 Aug 2025 21:05:30 -0400 Subject: [PATCH 12/13] comment out output tree --- offline/packages/tpc/TpcCombinedRawDataUnpacker.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/offline/packages/tpc/TpcCombinedRawDataUnpacker.h b/offline/packages/tpc/TpcCombinedRawDataUnpacker.h index 8eb269150a..9add5c69ae 100644 --- a/offline/packages/tpc/TpcCombinedRawDataUnpacker.h +++ b/offline/packages/tpc/TpcCombinedRawDataUnpacker.h @@ -144,7 +144,7 @@ class TpcCombinedRawDataUnpacker : public SubsysReco uint16_t t_sampch=std::numeric_limits::quiet_NaN(); - bool m_writeTree{true}; + bool m_writeTree{false}; bool m_do_baseline_corr{false}; int m_baseline_nsigma{2}; bool m_do_zs_emulation{false}; From 9ef42df57bdea10da4033db962e23edb93d48d48 Mon Sep 17 00:00:00 2001 From: Chenxi Ma <997873828mcx@gmail.com> Date: Wed, 3 Sep 2025 21:52:11 -0400 Subject: [PATCH 13/13] now it applies to both sides and include all pads within sigmas from center of the cloud --- .../g4tpc/PHG4TpcPadPlaneReadout.cc | 470 ++++++++++++------ .../g4tpc/PHG4TpcPadPlaneReadout.h | 18 +- 2 files changed, 342 insertions(+), 146 deletions(-) diff --git a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc index ce78ba0b52..96bdff938d 100644 --- a/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc +++ b/simulation/g4simulation/g4tpc/PHG4TpcPadPlaneReadout.cc @@ -67,7 +67,7 @@ namespace return std::exp(-square(x / sigma) / 2) / (sigma * std::sqrt(2 * M_PI)); } - static constexpr unsigned int print_layer = 18; + static constexpr int print_layer = 18; } // namespace @@ -91,6 +91,43 @@ PHG4TpcPadPlaneReadout::~PHG4TpcPadPlaneReadout() { 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; + + PHG4TpcCylinderGeomContainer::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) +PHG4TpcCylinderGeom* +PHG4TpcPadPlaneReadout::getGeomForLayer(unsigned int layer) const +{ + PHG4TpcCylinderGeomContainer::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; } //_________________________________________________________ @@ -402,14 +439,36 @@ double PHG4TpcPadPlaneReadout::getSingleEGEMAmplification(TF1 *f) } -//_________________________________________________________ +/* //_________________________________________________________ inline void rotatePointToSector(double x, double y, double& xNew, double& yNew, - int& sector) + int& sector, + const unsigned int side) { // 1) compute original phi in [−π, +π] double phi = std::atan2(y, x); + sector = -1; + for (int s = 0; s < 12; ++s) { + double min_phi = sector_min_Phi[side][s]; + double max_phi = sector_max_Phi[side][s]; + + // Check if phi is in this sector + // Note: need to handle wraparound at ±π + if (min_phi <= max_phi) { + if (phi >= min_phi && phi <= max_phi) { + sector = s; + break; + } + } else { // wraps around ±π + if (phi >= min_phi || phi <= max_phi) { + sector = s; + break; + } + } + } + + double dphi = sector_min_phi[side][sector] - sector_min_phi[side][2]; // 2) find the 30°‐wide wedge it lives in const double PI = std::acos(-1.0); double wedgeWidth = 2.0 * PI / TpcDefs::NSectors; // = π/6 @@ -429,10 +488,86 @@ inline void rotatePointToSector(double x, double y, double phiRot = phi + dphi; xNew = R * std::cos(phiRot); yNew = R * std::sin(phiRot); +} */ +void PHG4TpcPadPlaneReadout::rotatePointToSector( + double x, double y, + unsigned int side, + int& sectorFound, + double& xNew, double& yNew +) +{ + // ---- 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 sector center to 12 o'clock ------------------------------- + const double centerPhi = mid(sector_min_Phi[side][sectorFound], + sector_max_Phi[side][sectorFound]); + const double dphi = wrap(PI/2.0 - centerPhi); + const double phiRot = wrap(phi + dphi); + + xNew = R * std::cos(phiRot); + yNew = R * std::sin(phiRot); + + // ---- 3) mirror for North side to reuse South polygons -------------------- + + constexpr unsigned NORTH_SIDE = 1; // + if (side == NORTH_SIDE) + { + // reflect left/right in the canonical frame + xNew = -xNew; + // yNew unchanged + } } + /* //_________________________________________________________ std::vector computeShaperKernel() { @@ -590,10 +725,8 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( } // store phi bins and tbins upfront to avoid repetitive checks on the phi methods - const auto phibins = LayerGeom->get_phibins(); - /* pass_data.nphibins = phibins; */ - - const auto tbins = LayerGeom->get_zbins(); + /* const auto phibins = LayerGeom->get_phibins(); */ + /* const auto tbins = LayerGeom->get_zbins(); */ sector_min_Phi = LayerGeom->get_sector_min_phi(); sector_max_Phi = LayerGeom->get_sector_max_phi(); @@ -612,6 +745,16 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( //=============================== double nelec = getSingleEGEMAmplification(); + + // consider all layers whose radial annulus intersects a 5σ window + const double nsig_window = 5.0; + + std::vector cand_layers = layersInRadialWindow(rad_gem, sigmaT, nsig_window); + if (cand_layers.empty()) { + cand_layers.push_back(layernum); + } + + // Applying weight with respect to the rad_gem and phi after electrons are redistributed double phi_gain = phi; if (phi < 0) @@ -710,37 +853,10 @@ void PHG4TpcPadPlaneReadout::MapToPadPlane( - std::vector pad_phibin; - std::vector pad_phibin_share; - - //std::vector pad_phibin_ref; - //std::vector pad_phibin_share_ref; - //std::vector pad_phibin_serf; - //std::vector pad_phibin_share_serf; - - // populate_zigzag_phibins(side, layernum, phi, sigmaT, pad_phibin_ref, pad_phibin_share_ref); - - - populate_zigzag_phibins(side, layernum, phi, sigmaT, pad_phibin, pad_phibin_share); - - //SERF_zigzag_phibins(side, layernum, phi, rad_gem, sigmaT, 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) - { - double pad_share = pad_phibin_share[ipad]; - norm1 += pad_share; - } - for (unsigned int iphi = 0; iphi < pad_phibin.size(); ++iphi) - { - pad_phibin_share[iphi] /= norm1; - } + // 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) @@ -789,109 +905,188 @@ norm1 = 0.0; { adc_tbin_share[it] /= tnorm; } - - // 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) + if (m_use_serf_padsharing) { - int pad_num = pad_phibin[ipad]; - double pad_share = pad_phibin_share[ipad]; - - for (unsigned int it = 0; it < adc_tbin.size(); ++it) + // Pass 1: compute total Gaussian mass across all candidate layers + double total_mass = 0.0; + for (unsigned int layer_cand : cand_layers) { - int tbin_num = adc_tbin[it]; - double adc_bin_share = adc_tbin_share[it]; + PHG4TpcCylinderGeom* 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; + } - // 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 (total_mass <= 1e-16) total_mass = 1.0; // avoid division by zero - if (tbin_num >= tbins) - { - std::cout << " Error making key: adc_tbin " << tbin_num << " ntbins " << tbins << std::endl; - } - if (pad_num >= phibins) + // Pass 2: fill hits with correct per-layer mass fraction + for (unsigned int layer_cand : cand_layers) + { + PHG4TpcCylinderGeom* 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 << " Error making key: pad_phibin " << pad_num << " nphibins " << phibins << 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; - // 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) - { - std::cout << " tbin_num " << tbin_num << " tcenter " << tcenter << " pad_num " << pad_num << " phicenter " << phicenter - << " neffelectrons " << neffelectrons << " neffelectrons_threshold " << neffelectrons_threshold << std::endl; - } + 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); - // 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 - - // The side is an input parameter - - // 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(layernum, 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); - - // generate the key for this hit, requires tbin and phibin - TrkrDefs::hitkey hitkey = TpcDefs::genHitKey((unsigned int) pad_num, (unsigned int) tbin_num); - // See if this hit already exists - TrkrHit *hit = nullptr; - hit = hitsetit->second->getHit(hitkey); - if (!hit) - { - // create a new one - hit = new TrkrHitv2(); - hitsetit->second->addHitSpecificKey(hitkey, hit); + 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); + + tpc_truth_clusterer.addhitset(hitsetkey, hitkey, neffelectrons_bin); + + 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; + } } - // Either way, add the energy to it -- adc values will be added at digitization - hit->addEnergy(neffelectrons); - - tpc_truth_clusterer.addhitset(hitsetkey, hitkey, neffelectrons); + } + } + else + { + // 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) + { + PHG4TpcCylinderGeom* 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; - // 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) + // Pass 2: per-layer phi sharing and fill + for (unsigned int layer_cand : cand_layers) + { + PHG4TpcCylinderGeom* 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 - 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); + 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; - /* - if (Verbosity() > 0) - { - assert(nthit); - nthit->Fill(layernum, pad_num, tbin_num, neffelectrons); - } - */ + 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); - } // end of loop over adc T bins - } // end of loop over zigzag pads + 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); + + tpc_truth_clusterer.addhitset(hitsetkey, hitkey, neffelectrons_bin); + + 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 pad sharing and hit filling across candidate layers /* pass_data.phi_integral = phi_integral; */ /* pass_data.time_integral = t_integral; */ @@ -904,22 +1099,11 @@ norm1 = 0.0; } */ - 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++; @@ -1078,7 +1262,7 @@ void PHG4TpcPadPlaneReadout::SERF_zigzag_phibins(const unsigned int side, const 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 = "<