File indexing completed on 2026-08-30 08:14:28
0001 #include "TpcCentralMembraneStripeMatching.h"
0002 #include "StripeComparison.h"
0003 #include "StripeDetector.h"
0004 #include "helpers.h"
0005 #include "parameters.h"
0006
0007 #include <cdbobjects/CDBTTree.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010
0011 #include <phool/getClass.h>
0012 #include <phool/PHCompositeNode.h>
0013
0014 #include <trackbase/LaserCluster.h>
0015 #include <trackbase/LaserClusterContainer.h>
0016 #include <trackbase/TpcDefs.h>
0017
0018 #include <TFile.h>
0019 #include <TDirectory.h>
0020 #include <TH2.h>
0021 #include <TH3.h>
0022 #include <TVector3.h>
0023
0024 #include <algorithm>
0025 #include <cmath>
0026 #include <iostream>
0027 #include <string>
0028 #include <utility>
0029 #include <vector>
0030
0031
0032 TpcCentralMembraneStripeMatching::TpcCentralMembraneStripeMatching(const std::string &name) : SubsysReco(name) {}
0033
0034 int TpcCentralMembraneStripeMatching::InitRun(PHCompositeNode * )
0035 {
0036 if (m_fillReferenceHistogramsOnly)
0037 {
0038 CreateReferenceHistograms();
0039 return 0;
0040 }
0041
0042 if (!useIdealStripesAsReference)
0043 {
0044 LoadReferenceHistograms();
0045 }
0046 CreateMeasuredHistograms();
0047 return 0;
0048 }
0049
0050 int TpcCentralMembraneStripeMatching::process_event(PHCompositeNode *topNode)
0051 {
0052 const int nodeStatus = GetNodes(topNode);
0053 if (nodeStatus != Fun4AllReturnCodes::EVENT_OK)
0054 {
0055 return nodeStatus;
0056 }
0057
0058 if (m_fillReferenceHistogramsOnly)
0059 {
0060 FillReferenceHistograms();
0061 }
0062 else
0063 {
0064 FillMeasuredHistograms();
0065 }
0066 return Fun4AllReturnCodes::EVENT_OK;
0067 }
0068
0069 int TpcCentralMembraneStripeMatching::End(PHCompositeNode * )
0070 {
0071 if (m_fillReferenceHistogramsOnly)
0072 {
0073 std::cout << "Filled reference hPetal histograms from LASER_CLUSTER over " << m_processedEvents << " events: negz=" << m_filledReferenceClusters[0] << " posz=" << m_filledReferenceClusters[1] << std::endl;
0074 return Fun4AllReturnCodes::EVENT_OK;
0075 }
0076
0077 TFile *outputfile = new TFile(m_outputfile.c_str(), "RECREATE");
0078 std::cout << "Writing output to file: " << m_outputfile << std::endl;
0079 std::cout << "Filled measured hPetal histograms from LASER_CLUSTER over " << m_processedEvents << " events: negz=" << m_filledClusters[0] << " posz=" << m_filledClusters[1] << std::endl;
0080 outputfile->cd();
0081
0082 if (useIdealStripesAsReference)
0083 {
0084
0085 CDBTTree *cdbttree = new CDBTTree(m_idealStripePatternFile);
0086 cdbttree->LoadCalibrations();
0087 auto cdbMap = cdbttree->GetDoubleEntryMap();
0088 m_stripes[1][0].clear();
0089 m_stripes[1][1].clear();
0090 for (const auto &entry : cdbMap)
0091 {
0092 const auto index = entry.first;
0093 std::array<double, 3> stripe{};
0094 stripe[stripe_r] = cdbttree->GetDoubleValue(index, "truthR");
0095 stripe[stripe_phi] = cdbttree->GetDoubleValue(index, "truthPhi");
0096 if (stripe[stripe_phi] < 0.0)
0097 {
0098 stripe[stripe_phi] += 2 * M_PI;
0099 }
0100
0101 if (index > 180000)
0102 {
0103 m_stripes[1][0].push_back(stripe);
0104 }
0105 else
0106 {
0107 m_stripes[1][1].push_back(stripe);
0108 }
0109 }
0110 delete cdbttree;
0111
0112 if (m_stripes[1][0].empty() || m_stripes[1][1].empty())
0113 {
0114 std::cerr << "Ideal stripe pattern file has no stripes on one side: " << m_idealStripePatternFile << std::endl;
0115 return Fun4AllReturnCodes::ABORTRUN;
0116 }
0117 std::cout << "Loaded ideal stripe positions from " << m_idealStripePatternFile << ": negz=" << m_stripes[1][0].size() << " posz=" << m_stripes[1][1].size() << std::endl;
0118 }
0119
0120
0121
0122 if (CleanClusterHistograms() != Fun4AllReturnCodes::EVENT_OK)
0123 {
0124 std::cout << "Error in CleanClusterHistograms" << std::endl;
0125 return Fun4AllReturnCodes::ABORTEVENT;
0126 }
0127
0128 if (!useIdealStripesAsReference)
0129 {
0130 for (int side = 0; side < 2; ++side)
0131 {
0132 m_stripes[1][side].clear();
0133 m_detector.detect(m_hPetalCleaned[1][side], m_stripes[1][side]);
0134 }
0135 }
0136
0137 for (int side = 0; side < 2; ++side)
0138 {
0139 m_stripes[0][side].clear();
0140 m_detector.detect(m_hPetalCleaned[0][side], m_stripes[0][side]);
0141 ComputeStripeComparisonMaps(m_stripes[0][side], m_stripes[1][side], m_lamination_radial_gap_centers[0][side], side);
0142 }
0143
0144 outputfile->Write();
0145 outputfile->Close();
0146
0147 return Fun4AllReturnCodes::EVENT_OK;
0148 }
0149
0150 void TpcCentralMembraneStripeMatching::setOutputfile(const std::string &outputfile)
0151 {
0152 m_outputfile = outputfile;
0153 }
0154
0155 TH2 *TpcCentralMembraneStripeMatching::cloneReferenceHistogram(int side, const std::string &name) const
0156 {
0157 if (side < 0 || side >= 2 || !m_hPetal[1][side])
0158 {
0159 return nullptr;
0160 }
0161
0162 std::string cloneName = name;
0163 if (cloneName.empty())
0164 {
0165 cloneName = std::string(m_hPetal[1][side]->GetName()) + "_clone";
0166 }
0167
0168 auto *clone = dynamic_cast<TH2 *>(m_hPetal[1][side]->Clone(cloneName.c_str()));
0169 if (clone)
0170 {
0171 clone->SetDirectory(nullptr);
0172 }
0173 return clone;
0174 }
0175
0176 int TpcCentralMembraneStripeMatching::GetNodes(PHCompositeNode *topNode)
0177 {
0178 m_laserClusterContainer = findNode::getClass<LaserClusterContainer>(topNode, "LASER_CLUSTER");
0179 if (!m_laserClusterContainer)
0180 {
0181 std::cout << "TpcCentralMembraneStripeMatching::GetNodes - LASER_CLUSTER node missing" << std::endl;
0182 return Fun4AllReturnCodes::ABORTRUN;
0183 }
0184
0185 m_dcc_in_module_edge = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerModuleEdge");
0186
0187
0188
0189
0190
0191 m_dcc_in_static = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerStatic");
0192
0193
0194
0195
0196
0197 return Fun4AllReturnCodes::EVENT_OK;
0198 }
0199
0200 void TpcCentralMembraneStripeMatching::CreateMeasuredHistograms()
0201 {
0202 const char *names[2] = {"hPetal_measured_negz_raw", "hPetal_measured_posz_raw"};
0203
0204 for (int side = 0; side < 2; ++side)
0205 {
0206 delete m_hPetal[0][side];
0207 m_hPetal[0][side] = new TH2D(names[side], names[side], m_phiBins, 0.0, 2.0 * M_PI, m_rBins, m_rMin, m_rMax);
0208 m_hPetal[0][side]->SetDirectory(nullptr);
0209 }
0210 }
0211
0212 void TpcCentralMembraneStripeMatching::CreateReferenceHistograms()
0213 {
0214 const char *names[2] = {"hPetal_reference_negz_raw", "hPetal_reference_posz_raw"};
0215
0216 for (int side = 0; side < 2; ++side)
0217 {
0218 delete m_hPetal[1][side];
0219 m_hPetal[1][side] = new TH2D(names[side], names[side], m_phiBins, 0.0, 2.0 * M_PI, m_rBins, m_rMin, m_rMax);
0220 m_hPetal[1][side]->SetDirectory(nullptr);
0221 }
0222 }
0223
0224 void TpcCentralMembraneStripeMatching::LoadReferenceHistograms()
0225 {
0226 auto *referenceFile = TFile::Open(m_referenceStripePatternFile.c_str(), "READ");
0227 if (!referenceFile || referenceFile->IsZombie())
0228 {
0229 std::cout << "ERROR: Could not open reference stripe pattern file: " << m_referenceStripePatternFile << std::endl;
0230 return;
0231 }
0232 m_hPetal[1][0] = load_detached_histogram(referenceFile, "hPetal_South");
0233 m_hPetal[1][1] = load_detached_histogram(referenceFile, "hPetal_North");
0234 referenceFile->Close();
0235
0236 if (!m_hPetal[1][0] || !m_hPetal[1][1])
0237 {
0238 std::cout << "ERROR: Missing hPetal_South or hPetal_North reference histogram" << std::endl;
0239 }
0240 }
0241
0242 void TpcCentralMembraneStripeMatching::FillClusterHistograms(int histogramSet)
0243 {
0244 if (histogramSet < 0 || histogramSet >= 2 || !m_laserClusterContainer || !m_hPetal[histogramSet][0] || !m_hPetal[histogramSet][1])
0245 {
0246 return;
0247 }
0248
0249 ++m_processedEvents;
0250 const auto clusterRange = m_laserClusterContainer->getClusters();
0251 for (auto clusterIter = clusterRange.first; clusterIter != clusterRange.second; ++clusterIter)
0252 {
0253 const auto &[clusterKey, cluster] = *clusterIter;
0254 if (!cluster)
0255 {
0256 continue;
0257 }
0258
0259 bool side = (bool) TpcDefs::getSide(clusterKey);
0260 Acts::Vector3 pos(cluster->getX(), cluster->getY(), (side ? 1.0 : -1.0));
0261 if (m_dcc_in_module_edge)
0262 {
0263 pos = m_distortionCorrection.get_corrected_position(pos, m_dcc_in_module_edge);
0264 }
0265
0266 TVector3 tmp_pos(pos[0], pos[1], pos[2]);
0267
0268
0269 const unsigned int nLayers = cluster->getNLayers();
0270 if (nLayers < 2 || cluster->getSDWeightedLayer() >= 0.5)
0271 {
0272 continue;
0273 }
0274
0275 const int iside = side ? 1 : 0;
0276
0277 m_hPetal[histogramSet][iside]->Fill(NormalizeClusterPhi(tmp_pos.Phi()), tmp_pos.Perp());
0278
0279 if (histogramSet == 0)
0280 {
0281 ++m_filledClusters[iside];
0282 }
0283 else
0284 {
0285 ++m_filledReferenceClusters[iside];
0286 }
0287 }
0288 }
0289
0290 void TpcCentralMembraneStripeMatching::FillMeasuredHistograms()
0291 {
0292 FillClusterHistograms(0);
0293 }
0294
0295 void TpcCentralMembraneStripeMatching::FillReferenceHistograms()
0296 {
0297 FillClusterHistograms(1);
0298 }
0299
0300 double TpcCentralMembraneStripeMatching::NormalizeClusterPhi(double phi) const
0301 {
0302 while (phi < 0.0)
0303 {
0304 phi += 2.0 * M_PI;
0305 }
0306 while (phi >= 2.0 * M_PI)
0307 {
0308 phi -= 2.0 * M_PI;
0309 }
0310 return phi;
0311 }
0312
0313 void TpcCentralMembraneStripeMatching::ComputeStripeComparisonMaps(const std::vector<std::array<double, 3>> &measured, const std::vector<std::array<double, 3>> &reference, const std::vector<double> &measuredRadialGapCenters, int side)
0314 {
0315 if (!m_comparison.initialize(measured, reference, side, measuredRadialGapCenters) || !m_comparison.filter_isolated_inputs(measured, reference) || !m_comparison.build_global_pattern_matches() || !m_comparison.build_global_field_estimates())
0316 {
0317 m_comparison.clear();
0318 return;
0319 }
0320
0321 m_comparison.write_output_maps();
0322 m_comparison.write_corrected_measured_histogram(m_hPetal[0][side]);
0323 WriteStaticCorrectedMeasuredHistogram(side);
0324 m_comparison.write_distorted_reference_histogram(m_hPetal[1][side]);
0325 m_comparison.clear();
0326 }
0327
0328 void TpcCentralMembraneStripeMatching::WriteStaticCorrectedMeasuredHistogram(int side)
0329 {
0330 if (side < 0 || side >= 2 || !m_hPetal[0][side])
0331 {
0332 return;
0333 }
0334
0335 TDirectory *outputDirectory = gDirectory;
0336 const char *sideNames[2] = {"negz", "posz"};
0337 auto *staticFile = TFile::Open(m_staticCorrectionFile.c_str(), "READ");
0338 if (!staticFile || staticFile->IsZombie())
0339 {
0340 std::cout << "WARNING: Could not open static correction map file: " << m_staticCorrectionFile << std::endl;
0341 if (outputDirectory)
0342 {
0343 outputDirectory->cd();
0344 }
0345 return;
0346 }
0347
0348 const std::string suffix = sideNames[side];
0349 auto *deltaRMap = dynamic_cast<TH3 *>(staticFile->Get((std::string("hIntDistortionR_") + suffix).c_str()));
0350 auto *deltaPhiMap = dynamic_cast<TH3 *>(staticFile->Get((std::string("hIntDistortionP_") + suffix).c_str()));
0351 if (!deltaRMap || !deltaPhiMap)
0352 {
0353 std::cout << "WARNING: Missing hIntDistortionR_" << suffix << " or hIntDistortionP_" << suffix << " in static correction map file: " << m_staticCorrectionFile << std::endl;
0354 staticFile->Close();
0355 if (outputDirectory)
0356 {
0357 outputDirectory->cd();
0358 }
0359 return;
0360 }
0361
0362 const auto closestZBinToCentralMembrane = [](TH3 *histogram) {
0363 int closestBin = 1;
0364 double closestAbsZ = std::abs(histogram->GetZaxis()->GetBinCenter(closestBin));
0365 for (int zBin = 2; zBin <= histogram->GetNbinsZ(); ++zBin)
0366 {
0367 const double absZ = std::abs(histogram->GetZaxis()->GetBinCenter(zBin));
0368 if (absZ < closestAbsZ)
0369 {
0370 closestAbsZ = absZ;
0371 closestBin = zBin;
0372 }
0373 }
0374 return closestBin;
0375 };
0376 const int zBinR = closestZBinToCentralMembrane(deltaRMap);
0377 const int zBinPhi = closestZBinToCentralMembrane(deltaPhiMap);
0378 std::cout << "Writing static corrected measured histogram for " << suffix
0379 << " using static map z bins R=" << zBinR << " (z=" << deltaRMap->GetZaxis()->GetBinCenter(zBinR)
0380 << ") P=" << zBinPhi << " (z=" << deltaPhiMap->GetZaxis()->GetBinCenter(zBinPhi) << ", interpreted as R#Delta#phi)" << std::endl;
0381
0382 TH2 *sourceHistogram = m_hPetal[0][side];
0383 const int nPhiBins = sourceHistogram->GetNbinsX();
0384 const int nRBins = sourceHistogram->GetNbinsY();
0385 auto *shifted = new TH2D((std::string("hPetal_measured_corrected_static_") + suffix).c_str(),
0386 (std::string("Measured cluster histogram shifted by static distortion map z=0 slice - ") + suffix + ";#phi [rad];R [cm]").c_str(),
0387 nPhiBins, sourceHistogram->GetXaxis()->GetXmin(), sourceHistogram->GetXaxis()->GetXmax(),
0388 nRBins, sourceHistogram->GetYaxis()->GetXmin(), sourceHistogram->GetYaxis()->GetXmax());
0389 shifted->SetDirectory(nullptr);
0390
0391 const double phiMin = sourceHistogram->GetXaxis()->GetXmin();
0392 const double phiMax = sourceHistogram->GetXaxis()->GetXmax();
0393 const double phiWidth = phiMax - phiMin;
0394 for (int phiBin = 1; phiBin <= nPhiBins; ++phiBin)
0395 {
0396 const double phi = sourceHistogram->GetXaxis()->GetBinCenter(phiBin);
0397 for (int rBin = 1; rBin <= nRBins; ++rBin)
0398 {
0399 const double content = sourceHistogram->GetBinContent(phiBin, rBin);
0400 if (content == 0.0)
0401 {
0402 continue;
0403 }
0404
0405 const double r = sourceHistogram->GetYaxis()->GetBinCenter(rBin);
0406 const int staticPhiBinR = deltaRMap->GetXaxis()->FindBin(phi);
0407 const int staticRBinR = deltaRMap->GetYaxis()->FindBin(r);
0408 const int staticPhiBinP = deltaPhiMap->GetXaxis()->FindBin(phi);
0409 const int staticRBinP = deltaPhiMap->GetYaxis()->FindBin(r);
0410 const double deltaR = deltaRMap->GetBinContent(staticPhiBinR, staticRBinR, zBinR);
0411 const double rDeltaPhi = deltaPhiMap->GetBinContent(staticPhiBinP, staticRBinP, zBinPhi);
0412 const double deltaPhi = std::abs(r) > 1e-6 ? rDeltaPhi / r : 0.0;
0413
0414 const double shiftedR = r - deltaR;
0415 double shiftedPhi = phi - deltaPhi;
0416 while (shiftedPhi < phiMin)
0417 {
0418 shiftedPhi += phiWidth;
0419 }
0420 while (shiftedPhi >= phiMax)
0421 {
0422 shiftedPhi -= phiWidth;
0423 }
0424
0425 const int shiftedPhiBin = shifted->GetXaxis()->FindBin(shiftedPhi);
0426 const int shiftedRBin = shifted->GetYaxis()->FindBin(shiftedR);
0427 if (shiftedPhiBin < 1 || shiftedPhiBin > nPhiBins || shiftedRBin < 1 || shiftedRBin > nRBins)
0428 {
0429 continue;
0430 }
0431 shifted->SetBinContent(shiftedPhiBin, shiftedRBin, shifted->GetBinContent(shiftedPhiBin, shiftedRBin) + content);
0432 }
0433 }
0434
0435 shifted->SetStats(false);
0436 if (outputDirectory)
0437 {
0438 outputDirectory->cd();
0439 }
0440 safe_write_object(shifted);
0441 delete shifted;
0442 staticFile->Close();
0443 if (outputDirectory)
0444 {
0445 outputDirectory->cd();
0446 }
0447 }
0448
0449 int TpcCentralMembraneStripeMatching::CleanClusterHistograms()
0450 {
0451 if (!m_hPetal[0][0] || !m_hPetal[0][1])
0452 {
0453 std::cout << "WARNING: No measured stripe pattern histograms filled" << std::endl;
0454 return Fun4AllReturnCodes::ABORTEVENT;
0455 }
0456
0457 if (!useIdealStripesAsReference && (!m_hPetal[1][0] || !m_hPetal[1][1]))
0458 {
0459 std::cout << "WARNING: No reference stripe pattern histograms loaded" << std::endl;
0460 return Fun4AllReturnCodes::ABORTEVENT;
0461 }
0462
0463 const char *run_names[2] = {"measured", "reference"};
0464 const char *side_names[2] = {"negz", "posz"};
0465 const int run_count = useIdealStripesAsReference ? 1 : 2;
0466 for (int r = 0; r < run_count; r++)
0467 {
0468 for (int s = 0; s < 2; s++)
0469 {
0470 if (m_hPetal[r][s])
0471 {
0472 m_hPetal[r][s]->SetDirectory(nullptr);
0473 m_hPetal[r][s]->SetName((std::string("hPetal_") + run_names[r] + "_" + side_names[s]).c_str());
0474 safe_write_object(m_hPetal[r][s]);
0475 }
0476
0477 if (!applyPreStripeLaminationMask)
0478 {
0479 delete m_hPetalCleaned[r][s];
0480 m_hPetalCleaned[r][s] = dynamic_cast<TH2 *>(m_hPetal[r][s]->Clone((std::string(m_hPetal[r][s]->GetName()) + "_cleaned").c_str()));
0481 if (m_hPetalCleaned[r][s])
0482 {
0483 m_hPetalCleaned[r][s]->SetDirectory(nullptr);
0484 safe_write_object(m_hPetalCleaned[r][s]);
0485 }
0486 else
0487 {
0488 std::cout << "ERROR: Failed to clone unmasked histogram for " << run_names[r] << " " << side_names[s] << std::endl;
0489 return Fun4AllReturnCodes::ABORTEVENT;
0490 }
0491 m_lamination_radial_gap_centers[r][s].clear();
0492 continue;
0493 }
0494
0495 delete m_hPetalCleaned[r][s];
0496 auto cleaning = clean_laminations(m_hPetal[r][s]);
0497 m_hPetalCleaned[r][s] = cleaning.first;
0498 m_lamination_radial_gap_centers[r][s] = std::move(cleaning.second);
0499
0500 if (!m_hPetalCleaned[r][s])
0501 {
0502 std::cout << "ERROR: Failed to build cleaned histogram for " << run_names[r] << " " << side_names[s] << std::endl;
0503 return Fun4AllReturnCodes::ABORTEVENT;
0504 }
0505 safe_write_object(m_hPetalCleaned[r][s]);
0506 }
0507 }
0508
0509 return Fun4AllReturnCodes::EVENT_OK;
0510 }
