Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:17:33

0001 #include "ClusterCDFCalculator.h"
0002 
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TMatrixD.h>
0006 #include <TVectorD.h>
0007 
0008 #include <boost/format.hpp>
0009 
0010 #include <iostream>
0011 
0012 ClusterCDFCalculator::~ClusterCDFCalculator()
0013 {
0014   if (file)
0015   {
0016     file->Close();
0017     delete file;
0018   }
0019 }
0020 
0021 void ClusterCDFCalculator::LoadProfile(const std::string &filename)
0022 {
0023   if (file)
0024   {
0025     std::cout << "ClusterCDFCalculator::LoadProfile Warning: a file is already loaded.. closing it first" << std::endl;
0026     file->Close();
0027     delete file;
0028     file = nullptr;
0029   }
0030 
0031   file = new TFile(filename.c_str(), "read");
0032 
0033   if (!file || file->IsZombie())
0034   {
0035     std::cout << "ClusterCDFCalculator::LoadProfile error: Could not open file!" << std::endl;
0036     return;
0037   }
0038 
0039   LoadHistogramsAndMatrices();
0040 }
0041 
0042 void ClusterCDFCalculator::LoadHistogramsAndMatrices()
0043 {
0044   file->GetObject("h_photon_hen", henbins_photon);
0045   if (!henbins_photon)
0046   {
0047     std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() FATAL ERROR: photon energy bins histogram does not exist.. return.." << std::endl;
0048     return;
0049   }
0050 
0051   file->GetObject("h_pi0_hen", henbins_pi0);
0052   if (!henbins_pi0)
0053   {
0054     std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() FATAL ERROR: pi0 energy bins histogram does not exist.. return.." << std::endl;
0055     return;
0056   }
0057 
0058   if (henbins_photon->GetNbinsX() != henbins_photon->GetNbinsX())
0059   {
0060     std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() FATAL ERROR: inconsistent energy bins for photon and pi0... return.." << std::endl;
0061     return;
0062   }
0063 
0064   const int nBins = henbins_photon->GetNbinsX();
0065 
0066   // Resize photon
0067   hD2_3x3_photon.resize(nBins, nullptr);
0068   hD2_5x5_photon.resize(nBins, nullptr);
0069   hD2_7x7_photon.resize(nBins, nullptr);
0070   meanD2_3x3_photon.resize(nBins, -1);
0071   meanD2_5x5_photon.resize(nBins, -1);
0072   meanD2_7x7_photon.resize(nBins, -1);
0073   inverseCovarianceMatrix_3x3_photon.resize(nBins, TMatrixD(NMATRIX_3x3, NMATRIX_3x3));
0074   inverseCovarianceMatrix_5x5_photon.resize(nBins, TMatrixD(NMATRIX_5x5, NMATRIX_5x5));
0075   inverseCovarianceMatrix_7x7_photon.resize(nBins, TMatrixD(NMATRIX_7x7, NMATRIX_7x7));
0076   ratioHistograms_3x3_photon.resize(nBins, std::vector<TH1 *>(NMATRIX_3x3, nullptr));
0077   ratioHistograms_5x5_photon.resize(nBins, std::vector<TH1 *>(NMATRIX_5x5, nullptr));
0078   ratioHistograms_7x7_photon.resize(nBins, std::vector<TH1 *>(NMATRIX_7x7, nullptr));
0079 
0080   // Resize pi0
0081   hD2_3x3_pi0.resize(nBins, nullptr);
0082   hD2_5x5_pi0.resize(nBins, nullptr);
0083   hD2_7x7_pi0.resize(nBins, nullptr);
0084   meanD2_3x3_pi0.resize(nBins, -1);
0085   meanD2_5x5_pi0.resize(nBins, -1);
0086   meanD2_7x7_pi0.resize(nBins, -1);
0087   inverseCovarianceMatrix_3x3_pi0.resize(nBins, TMatrixD(NMATRIX_3x3, NMATRIX_3x3));
0088   inverseCovarianceMatrix_5x5_pi0.resize(nBins, TMatrixD(NMATRIX_5x5, NMATRIX_5x5));
0089   inverseCovarianceMatrix_7x7_pi0.resize(nBins, TMatrixD(NMATRIX_7x7, NMATRIX_7x7));
0090   ratioHistograms_3x3_pi0.resize(nBins, std::vector<TH1 *>(NMATRIX_3x3, nullptr));
0091   ratioHistograms_5x5_pi0.resize(nBins, std::vector<TH1 *>(NMATRIX_5x5, nullptr));
0092   ratioHistograms_7x7_pi0.resize(nBins, std::vector<TH1 *>(NMATRIX_7x7, nullptr));
0093 
0094   for (int binidx = 0; binidx < nBins; ++binidx)
0095   {
0096     // photon hist
0097     std::string hD2_3x3_name_photon = (boost::format("h_photon_hD2_3x3_en%d") % binidx).str();
0098     std::string hD2_5x5_name_photon = (boost::format("h_photon_hD2_5x5_en%d") % binidx).str();
0099     std::string hD2_7x7_name_photon = (boost::format("h_photon_hD2_7x7_en%d") % binidx).str();
0100 
0101     file->GetObject(hD2_3x3_name_photon.c_str(), hD2_3x3_photon[binidx]);
0102     file->GetObject(hD2_5x5_name_photon.c_str(), hD2_5x5_photon[binidx]);
0103     file->GetObject(hD2_7x7_name_photon.c_str(), hD2_7x7_photon[binidx]);
0104 
0105     if (!hD2_3x3_photon[binidx] || !hD2_5x5_photon[binidx] || !hD2_7x7_photon[binidx])
0106     {
0107       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: One or more photon hD2 histograms for bin " << binidx << " are missing. Returning..." << std::endl;
0108       return;
0109     }
0110 
0111     TH1 *hD2mean3_photon{nullptr};
0112     file->GetObject((boost::format("h_photon_hD2mean3_en%d") % binidx).str().c_str(), hD2mean3_photon);
0113     TH1 *hD2mean5_photon{nullptr};
0114     file->GetObject((boost::format("h_photon_hD2mean5_en%d") % binidx).str().c_str(), hD2mean5_photon);
0115     TH1 *hD2mean7_photon{nullptr};
0116     file->GetObject((boost::format("h_photon_hD2mean7_en%d") % binidx).str().c_str(), hD2mean7_photon);
0117 
0118     if (!hD2mean3_photon || !hD2mean5_photon || !hD2mean7_photon)
0119     {
0120       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: One or more required photon mean histograms are missing for bin " << binidx << std::endl;
0121       return;
0122     }
0123     meanD2_3x3_photon[binidx] = hD2mean3_photon->GetBinContent(1);
0124     meanD2_5x5_photon[binidx] = hD2mean5_photon->GetBinContent(1);
0125     meanD2_7x7_photon[binidx] = hD2mean7_photon->GetBinContent(1);
0126 
0127     for (int i = 0; i < NMATRIX_3x3; ++i)
0128     {
0129       std::string histName = (boost::format("h_photon_heratio_3x3_en%d_%d") % binidx % i).str();
0130       file->GetObject(histName.c_str(), ratioHistograms_3x3_photon[binidx][i]);
0131       if (!ratioHistograms_3x3_photon[binidx][i])
0132       {
0133         std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: photon hist " << histName.c_str() << " is missing." << std::endl;
0134         return;
0135       }
0136     }
0137     for (int i = 0; i < NMATRIX_5x5; ++i)
0138     {
0139       std::string histName = (boost::format("h_photon_heratio_5x5_en%d_%d") % binidx % i).str();
0140       file->GetObject(histName.c_str(), ratioHistograms_5x5_photon[binidx][i]);
0141       if (!ratioHistograms_5x5_photon[binidx][i])
0142       {
0143         std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: photon hist " << histName.c_str() << " is missing." << std::endl;
0144         return;
0145       }
0146     }
0147     for (int i = 0; i < NMATRIX_7x7; ++i)
0148     {
0149       std::string histName = (boost::format("h_photon_heratio_7x7_en%d_%d") % binidx % i).str();
0150       file->GetObject(histName.c_str(), ratioHistograms_7x7_photon[binidx][i]);
0151       if (!ratioHistograms_7x7_photon[binidx][i])
0152       {
0153         std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: photon hist " << histName.c_str() << " is missing." << std::endl;
0154         return;
0155       }
0156     }
0157     TH2 *hCovMatrix3x3_photon{nullptr};
0158     file->GetObject((boost::format("h_photon_hCovMatrix3_en%d") % binidx).str().c_str(), hCovMatrix3x3_photon);
0159     TH2 *hCovMatrix5x5_photon{nullptr};
0160     file->GetObject((boost::format("h_photon_hCovMatrix5_en%d") % binidx).str().c_str(), hCovMatrix5x5_photon);
0161     TH2 *hCovMatrix7x7_photon{nullptr};
0162     file->GetObject((boost::format("h_photon_hCovMatrix7_en%d") % binidx).str().c_str(), hCovMatrix7x7_photon);
0163 
0164     if (!hCovMatrix3x3_photon || !hCovMatrix5x5_photon || !hCovMatrix7x7_photon)
0165     {
0166       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: photon covariance matrix histograms are missing for bin " << binidx << std::endl;
0167       return;
0168     }
0169     for (int i = 0; i < NMATRIX_3x3; ++i)
0170     {
0171       for (int j = 0; j < NMATRIX_3x3; ++j)
0172       {
0173         inverseCovarianceMatrix_3x3_photon[binidx](i, j) = hCovMatrix3x3_photon->GetBinContent(i + 1, j + 1);
0174       }
0175     }
0176     for (int i = 0; i < NMATRIX_5x5; ++i)
0177     {
0178       for (int j = 0; j < NMATRIX_5x5; ++j)
0179       {
0180         inverseCovarianceMatrix_5x5_photon[binidx](i, j) = hCovMatrix5x5_photon->GetBinContent(i + 1, j + 1);
0181       }
0182     }
0183     for (int i = 0; i < NMATRIX_7x7; ++i)
0184     {
0185       for (int j = 0; j < NMATRIX_7x7; ++j)
0186       {
0187         inverseCovarianceMatrix_7x7_photon[binidx](i, j) = hCovMatrix7x7_photon->GetBinContent(i + 1, j + 1);
0188       }
0189     }
0190     try
0191     {
0192       inverseCovarianceMatrix_3x3_photon[binidx].Invert();
0193       inverseCovarianceMatrix_5x5_photon[binidx].Invert();
0194       inverseCovarianceMatrix_7x7_photon[binidx].Invert();
0195     }
0196     catch (const std::exception &e)
0197     {
0198       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices - Exception in photon covariance invert for " << e.what() << " abort.." << std::endl;
0199       return;
0200     }
0201 
0202     // pi0 hist
0203     std::string hD2_3x3_name_pi0 = (boost::format("h_pi0_hD2_3x3_en%d") % binidx).str();
0204     std::string hD2_5x5_name_pi0 = (boost::format("h_pi0_hD2_5x5_en%d") % binidx).str();
0205     std::string hD2_7x7_name_pi0 = (boost::format("h_pi0_hD2_7x7_en%d") % binidx).str();
0206 
0207     file->GetObject(hD2_3x3_name_pi0.c_str(), hD2_3x3_pi0[binidx]);
0208     file->GetObject(hD2_5x5_name_pi0.c_str(), hD2_5x5_pi0[binidx]);
0209     file->GetObject(hD2_7x7_name_pi0.c_str(), hD2_7x7_pi0[binidx]);
0210 
0211     if (!hD2_3x3_pi0[binidx] || !hD2_5x5_pi0[binidx] || !hD2_7x7_pi0[binidx])
0212     {
0213       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: One or more pi0 hD2 histograms for bin " << binidx << " are missing. Returning..." << std::endl;
0214       return;
0215     }
0216 
0217     TH1 *hD2mean3_pi0{nullptr};
0218     file->GetObject((boost::format("h_pi0_hD2mean3_en%d") % binidx).str().c_str(), hD2mean3_pi0);
0219     TH1 *hD2mean5_pi0{nullptr};
0220     file->GetObject((boost::format("h_pi0_hD2mean5_en%d") % binidx).str().c_str(), hD2mean5_pi0);
0221     TH1 *hD2mean7_pi0{nullptr};
0222     file->GetObject((boost::format("h_pi0_hD2mean7_en%d") % binidx).str().c_str(), hD2mean7_pi0);
0223 
0224     if (!hD2mean3_pi0 || !hD2mean5_pi0 || !hD2mean7_pi0)
0225     {
0226       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: One or more required pi0 mean histograms are missing for bin " << binidx << std::endl;
0227       return;
0228     }
0229     meanD2_3x3_pi0[binidx] = hD2mean3_pi0->GetBinContent(1);
0230     meanD2_5x5_pi0[binidx] = hD2mean5_pi0->GetBinContent(1);
0231     meanD2_7x7_pi0[binidx] = hD2mean7_pi0->GetBinContent(1);
0232 
0233     for (int i = 0; i < NMATRIX_3x3; ++i)
0234     {
0235       std::string histName = (boost::format("h_pi0_heratio_3x3_en%d_%d") % binidx % i).str();
0236       file->GetObject(histName.c_str(), ratioHistograms_3x3_pi0[binidx][i]);
0237       if (!ratioHistograms_3x3_pi0[binidx][i])
0238       {
0239         std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: pi0 hist " << histName.c_str() << " is missing." << std::endl;
0240         return;
0241       }
0242     }
0243     for (int i = 0; i < NMATRIX_5x5; ++i)
0244     {
0245       std::string histName = (boost::format("h_pi0_heratio_5x5_en%d_%d") % binidx % i).str();
0246       file->GetObject(histName.c_str(), ratioHistograms_5x5_pi0[binidx][i]);
0247       if (!ratioHistograms_5x5_pi0[binidx][i])
0248       {
0249         std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: pi0 hist " << histName.c_str() << " is missing." << std::endl;
0250         return;
0251       }
0252     }
0253     for (int i = 0; i < NMATRIX_7x7; ++i)
0254     {
0255       std::string histName = (boost::format("h_pi0_heratio_7x7_en%d_%d") % binidx % i).str();
0256       file->GetObject(histName.c_str(), ratioHistograms_7x7_pi0[binidx][i]);
0257       if (!ratioHistograms_7x7_pi0[binidx][i])
0258       {
0259         std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: pi0 hist " << histName.c_str() << " is missing." << std::endl;
0260         return;
0261       }
0262     }
0263     TH2 *hCovMatrix3x3_pi0{nullptr};
0264     file->GetObject((boost::format("h_pi0_hCovMatrix3_en%d") % binidx).str().c_str(), hCovMatrix3x3_pi0);
0265     TH2 *hCovMatrix5x5_pi0{nullptr};
0266     file->GetObject((boost::format("h_pi0_hCovMatrix5_en%d") % binidx).str().c_str(), hCovMatrix5x5_pi0);
0267     TH2 *hCovMatrix7x7_pi0{nullptr};
0268     file->GetObject((boost::format("h_pi0_hCovMatrix7_en%d") % binidx).str().c_str(), hCovMatrix7x7_pi0);
0269 
0270     if (!hCovMatrix3x3_pi0 || !hCovMatrix5x5_pi0 || !hCovMatrix7x7_pi0)
0271     {
0272       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices() error: pi0 covariance matrix histograms are missing for bin " << binidx << std::endl;
0273       return;
0274     }
0275     for (int i = 0; i < NMATRIX_3x3; ++i)
0276     {
0277       for (int j = 0; j < NMATRIX_3x3; ++j)
0278       {
0279         inverseCovarianceMatrix_3x3_pi0[binidx](i, j) = hCovMatrix3x3_pi0->GetBinContent(i + 1, j + 1);
0280       }
0281     }
0282     for (int i = 0; i < NMATRIX_5x5; ++i)
0283     {
0284       for (int j = 0; j < NMATRIX_5x5; ++j)
0285       {
0286         inverseCovarianceMatrix_5x5_pi0[binidx](i, j) = hCovMatrix5x5_pi0->GetBinContent(i + 1, j + 1);
0287       }
0288     }
0289     for (int i = 0; i < NMATRIX_7x7; ++i)
0290     {
0291       for (int j = 0; j < NMATRIX_7x7; ++j)
0292       {
0293         inverseCovarianceMatrix_7x7_pi0[binidx](i, j) = hCovMatrix7x7_pi0->GetBinContent(i + 1, j + 1);
0294       }
0295     }
0296     try
0297     {
0298       inverseCovarianceMatrix_3x3_pi0[binidx].Invert();
0299       inverseCovarianceMatrix_5x5_pi0[binidx].Invert();
0300       inverseCovarianceMatrix_7x7_pi0[binidx].Invert();
0301     }
0302     catch (const std::exception &e)
0303     {
0304       std::cout << "ClusterCDFCalculator::LoadHistogramsAndMatrices - Exception in pi0 covariance invert for " << e.what() << " abort.." << std::endl;
0305       return;
0306     }
0307   }
0308 }
0309 
0310 double ClusterCDFCalculator::CalculateJointDeviation(const std::vector<double> &energies,
0311                                                      const std::vector<TH1 *> &ratioHistograms,
0312                                                      const TMatrixD &inverseCovarianceMatrix,
0313                                                      double meanD2)
0314 {
0315   int nDims = energies.size();
0316   double totalEnergy = 0.0;
0317   for (double energy : energies)
0318   {
0319     totalEnergy += energy;
0320   }
0321 
0322   TVectorD probVector(nDims);
0323   for (int i = 0; i < nDims; ++i)
0324   {
0325     double fraction = energies[i] / totalEnergy;
0326     int bin = std::min(std::max(ratioHistograms[i]->FindBin(fraction), 1), ratioHistograms[i]->GetNbinsX());
0327     double probability = std::max(ratioHistograms[i]->GetBinContent(bin), epsilon);
0328     probVector[i] = -log(probability);
0329   }
0330 
0331   return probVector * (inverseCovarianceMatrix * probVector) / meanD2;
0332 }
0333 
0334 double ClusterCDFCalculator::GetCDFValue(double jointDeviation, TH1 *hD2)
0335 {
0336   int bin = hD2->FindBin(jointDeviation);
0337   if (bin <= 0 || bin > hD2->GetNbinsX())
0338   {
0339     return 0.0;
0340   }
0341   return hD2->Integral(bin, hD2->GetNbinsX());
0342 }
0343 
0344 std::pair<double, double> ClusterCDFCalculator::GetCDF(const std::vector<double> &energies, double clusterenergy, int NMATRIXDIM)
0345 {
0346   if (NMATRIXDIM != 3 && NMATRIXDIM != 5 && NMATRIXDIM != 7)
0347   {
0348     std::cout << "ClusterCDFCalculator::GetCDF error: Invalid dimension as input! Only 3x3, 5x5, 7x7 are available (use input 3, 5, or 7)." << std::endl;
0349     return std::make_pair(-1, -1);
0350   }
0351 
0352   int towersize = energies.size();
0353   if (towersize != gridSize)
0354   {
0355     std::cout << "ClusterCDFCalculator::GetCDF error: Invalid tower energy vector size! It must be of size (49)." << std::endl;
0356     return std::make_pair(-1, -1);
0357   }
0358 
0359   int ibin = -1;
0360   if (henbins_photon)
0361   {
0362     ibin = henbins_photon->FindFixBin(clusterenergy) - 1;
0363     if (ibin < 0)
0364     {
0365       std::cout << "ClusterCDFCalculator::GetCDF fatal error: histogram bin below 0... set prob. both to -1." << std::endl;
0366       return std::make_pair(-1, -1);
0367     }
0368     if (ibin >= henbins_photon->GetNbinsX())
0369     {
0370       ibin = henbins_photon->GetNbinsX() - 1;
0371     }
0372   }
0373   else
0374   {
0375     std::cout << "ClusterCDFCalculator::GetCDF error: Photon energy bins histogram (henbins_photon) is not loaded!" << std::endl;
0376     return std::make_pair(-1, -1);
0377   }
0378 
0379   std::unordered_set<int> indices = getSubsetIndices(FULLGRIDSIZE, NMATRIXDIM);
0380   std::vector<double> inputenergies;
0381   for (int is = 0; is < towersize; ++is)
0382   {
0383     if (indices.count(is)) //NOLINT(readability-container-contains)
0384     {
0385       inputenergies.push_back(energies.at(is));
0386     }
0387   }
0388 
0389   double photonCDF = -1;
0390   double pi0CDF = -1;
0391   if (NMATRIXDIM == NMATRIXDIM3)
0392   {
0393     double jointDeviation_photon = CalculateJointDeviation(inputenergies, ratioHistograms_3x3_photon[ibin], inverseCovarianceMatrix_3x3_photon[ibin], meanD2_3x3_photon[ibin]);
0394     photonCDF = GetCDFValue(jointDeviation_photon, hD2_3x3_photon[ibin]);
0395 
0396     double jointDeviation_pi0 = CalculateJointDeviation(inputenergies, ratioHistograms_3x3_pi0[ibin], inverseCovarianceMatrix_3x3_pi0[ibin], meanD2_3x3_pi0[ibin]);
0397     pi0CDF = GetCDFValue(jointDeviation_pi0, hD2_3x3_pi0[ibin]);
0398   }
0399   else if (NMATRIXDIM == NMATRIXDIM5)
0400   {
0401     double jointDeviation_photon = CalculateJointDeviation(inputenergies, ratioHistograms_5x5_photon[ibin], inverseCovarianceMatrix_5x5_photon[ibin], meanD2_5x5_photon[ibin]);
0402     photonCDF = GetCDFValue(jointDeviation_photon, hD2_5x5_photon[ibin]);
0403 
0404     double jointDeviation_pi0 = CalculateJointDeviation(inputenergies, ratioHistograms_5x5_pi0[ibin], inverseCovarianceMatrix_5x5_pi0[ibin], meanD2_5x5_pi0[ibin]);
0405     pi0CDF = GetCDFValue(jointDeviation_pi0, hD2_5x5_pi0[ibin]);
0406   }
0407   else if (NMATRIXDIM == NMATRIXDIM7)
0408   {
0409     double jointDeviation_photon = CalculateJointDeviation(inputenergies, ratioHistograms_7x7_photon[ibin], inverseCovarianceMatrix_7x7_photon[ibin], meanD2_7x7_photon[ibin]);
0410     photonCDF = GetCDFValue(jointDeviation_photon, hD2_7x7_photon[ibin]);
0411 
0412     double jointDeviation_pi0 = CalculateJointDeviation(inputenergies, ratioHistograms_7x7_pi0[ibin], inverseCovarianceMatrix_7x7_pi0[ibin], meanD2_7x7_pi0[ibin]);
0413     pi0CDF = GetCDFValue(jointDeviation_pi0, hD2_7x7_pi0[ibin]);
0414   }
0415   else
0416   {
0417     std::cout << "ClusterCDFCalculator::GetCDF error: Invalid input size. Supported sizes are 3 (for 3x3), 5 (for 5x5), and 7 (for 7x7)." << std::endl;
0418     return std::make_pair(-1, -1);
0419   }
0420 
0421   return std::make_pair(photonCDF, pi0CDF);
0422 }
0423 
0424 std::unordered_set<int> ClusterCDFCalculator::getSubsetIndices(int _gridsize, int _subsetsize)
0425 {
0426   std::unordered_set<int> indices;
0427   if (_subsetsize > _gridsize || _subsetsize % 2 == 0)
0428   {
0429     std::cout << "Invalid subset size!" << std::endl;
0430     return indices;
0431   }
0432 
0433   int center = (_gridsize * _gridsize - 1) / 2;
0434   int halfSubset = _subsetsize / 2;
0435 
0436   for (int i = -halfSubset; i <= halfSubset; ++i)
0437   {
0438     for (int j = -halfSubset; j <= halfSubset; ++j)
0439     {
0440       int row = (center / _gridsize) + i;
0441       int col = (center % _gridsize) + j;
0442       if (row >= 0 && row < _gridsize && col >= 0 && col < _gridsize)
0443       {
0444         indices.insert((row * _gridsize) + col);
0445       }
0446     }
0447   }
0448   return indices;
0449 }