Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-12-16 09:20:28

0001 /**
0002  * \file TpcSpaceChargeReconstructionHelper.cc
0003  * \brief performs simple histogram manipulations for generating space charge distortion map suitable for correcting TPC clusters
0004  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0005  */
0006 
0007 #include "TpcSpaceChargeReconstructionHelper.h"
0008 
0009 #include <TH2.h>
0010 #include <TH3.h>
0011 #include <TString.h>
0012 
0013 #include <iostream>
0014 
0015 namespace
0016 {
0017   /// square
0018   template <class T>
0019   inline constexpr T square(T x)
0020   {
0021     return x * x;
0022   }
0023 
0024   // regularize angle between 0 and 2PI
0025   template <class T>
0026   inline constexpr T get_bound_angle(T phi)
0027   {
0028     while (phi < 0)
0029     {
0030       phi += 2 * M_PI;
0031     }
0032     while (phi >= 2 * M_PI)
0033     {
0034       phi -= 2 * M_PI;
0035     }
0036     return phi;
0037   }
0038 
0039   /// make sure angles in a given window are betwwen [0, 2PI]
0040   TpcSpaceChargeReconstructionHelper::range_t transform_range(const TpcSpaceChargeReconstructionHelper::range_t& a)
0041   {
0042     return {get_bound_angle(a.first), get_bound_angle(a.second)};
0043   }
0044 
0045   /// short class to check if a given value is in provided range
0046   class range_ftor_t
0047   {
0048    public:
0049     explicit range_ftor_t(double value)
0050       : m_value(value){};
0051 
0052     bool operator()(const TpcSpaceChargeReconstructionHelper::range_t& range) const
0053     {
0054       return m_value > range.first && m_value < range.second;
0055     }
0056 
0057    private:
0058     double m_value = 0;
0059   };
0060 
0061   /// returns true if given value is in provided range
0062   bool in_range(const double& value, const TpcSpaceChargeReconstructionHelper::range_t& range)
0063   { return range_ftor_t(value)(range); }
0064 
0065   /// returns true if given value is in any of the provided range
0066   bool in_range(const double& value, const TpcSpaceChargeReconstructionHelper::range_list_t& range_list)
0067   { return std::any_of(range_list.begin(), range_list.end(), range_ftor_t(value)); }
0068 
0069 }  // namespace
0070 
0071 /// default TPOT geometry, hardcoded
0072 /// phi range for central sector (TPC sectors 9 and 21)
0073 TpcSpaceChargeReconstructionHelper::range_t TpcSpaceChargeReconstructionHelper::phi_range_central = transform_range({-1.742, -1.43979});
0074 
0075 /// phi range for east sector (TPC sectors 8 and 20)
0076 TpcSpaceChargeReconstructionHelper::range_t TpcSpaceChargeReconstructionHelper::phi_range_east = transform_range({-2.27002, -1.9673});
0077 
0078 /// phi range for west sector (TPC sectors 10 and 22)
0079 TpcSpaceChargeReconstructionHelper::range_t TpcSpaceChargeReconstructionHelper::phi_range_west = transform_range({-1.21452, -0.911172});
0080 
0081 /// list of theta angles for each micromeas in central sector with theta defined as atan2(z,r)
0082 TpcSpaceChargeReconstructionHelper::range_list_t TpcSpaceChargeReconstructionHelper::theta_range_central = {{-0.918257, -0.613136}, {-0.567549, -0.031022}, {0.0332154, 0.570419}, {0.613631, 0.919122}};
0083 
0084 /// list of theta angles for each micromeas in central sector with theta defined as atan2(z,r)
0085 TpcSpaceChargeReconstructionHelper::range_list_t TpcSpaceChargeReconstructionHelper::theta_range_east = {{-0.636926, -0.133603}, {0.140678, 0.642714}};
0086 
0087 /// list of theta angles for each micromeas in central sector with theta defined as atan2(z,r)
0088 TpcSpaceChargeReconstructionHelper::range_list_t TpcSpaceChargeReconstructionHelper::theta_range_west = {{-0.643676, -0.141004}, {0.13485, 0.640695}};
0089 
0090 //____________________________________________________________________________________
0091 void TpcSpaceChargeReconstructionHelper::create_tpot_mask(TH3* hmask)
0092 {
0093   hmask->Reset();
0094 
0095   // loop over bins
0096   for (int ip = 0; ip < hmask->GetNbinsX(); ++ip)
0097   {
0098     for (int ir = 0; ir < hmask->GetNbinsY(); ++ir)
0099     {
0100       for (int iz = 0; iz < hmask->GetNbinsZ(); ++iz)
0101       {
0102         const double phi = hmask->GetXaxis()->GetBinCenter(ip + 1);
0103         const double r = hmask->GetYaxis()->GetBinCenter(ir + 1);
0104         const double z = hmask->GetZaxis()->GetBinCenter(iz + 1);
0105         const double theta = std::atan2(z, r);
0106         if (in_range(phi, phi_range_central) && in_range(theta, theta_range_central))
0107         {
0108           hmask->SetBinContent(ip + 1, ir + 1, iz + 1, 1);
0109         }
0110       }
0111     }
0112   }
0113 }
0114 
0115 //____________________________________________________________________________________
0116 void TpcSpaceChargeReconstructionHelper::extrapolate_z(TH3* source, const TH3* mask)
0117 {
0118   if (!source)
0119   {
0120     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_z - invalid source histogram" << std::endl;
0121     return;
0122   }
0123 
0124   if (!mask)
0125   {
0126     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_z - invalid mask histogram" << std::endl;
0127     return;
0128   }
0129 
0130   // loop over phi bins
0131   for (int ir = 0; ir < source->GetNbinsY(); ++ir)
0132   {
0133     for (int ip = 0; ip < source->GetNbinsX(); ++ip)
0134     {
0135       int iz_min = -1;
0136       for (int iz = 0; iz < source->GetNbinsZ(); ++iz)
0137       {
0138         bool in_range = mask->GetBinContent(ip + 1, ir + 1, iz + 1) > 0;
0139         if (in_range)
0140         {
0141           iz_min = iz;
0142           continue;
0143         }
0144 
0145         // iz is not in range. Check if iz_min was set. If not, skip this bin. If yes, find next bin in range
0146         if (iz_min < 0)
0147         {
0148           continue;
0149         }
0150 
0151         int iz_max = iz + 1;
0152         for (; iz_max < source->GetNbinsZ(); ++iz_max)
0153         {
0154           in_range = mask->GetBinContent(ip + 1, ir + 1, iz_max + 1) > 0;
0155           if (in_range)
0156           {
0157             break;
0158           }
0159         }
0160 
0161         // check interpolation range validity
0162         if (iz_max == source->GetNbinsZ())
0163         {
0164           continue;
0165         }
0166 
0167         // do the interpolation
0168         const double z_min = source->GetZaxis()->GetBinCenter(iz_min + 1);
0169         const double z_max = source->GetZaxis()->GetBinCenter(iz_max + 1);
0170         const double z = source->GetZaxis()->GetBinCenter(iz + 1);
0171         const double alpha_min = (z_max - z) / (z_max - z_min);
0172         const double alpha_max = (z - z_min) / (z_max - z_min);
0173 
0174         const double distortion_min = source->GetBinContent(ip + 1, ir + 1, iz_min + 1);
0175         const double distortion_max = source->GetBinContent(ip + 1, ir + 1, iz_max + 1);
0176         const double distortion = alpha_min * distortion_min + alpha_max * distortion_max;
0177 
0178         const double error_min = source->GetBinError(ip + 1, ir + 1, iz_min + 1);
0179         const double error_max = source->GetBinError(ip + 1, ir + 1, iz_max + 1);
0180         const double error = std::sqrt(square(alpha_min * error_min) + square(alpha_max * error_max));
0181 
0182         source->SetBinContent(ip + 1, ir + 1, iz + 1, distortion);
0183         source->SetBinError(ip + 1, ir + 1, iz + 1, error);
0184       }
0185     }
0186   }
0187 }
0188 
0189 //____________________________________________________________________________________
0190 void TpcSpaceChargeReconstructionHelper::extrapolate_z2(TH3* source, const TH3* mask, uint8_t side)
0191 {
0192   if (!source)
0193   {
0194     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_z - invalid source histogram" << std::endl;
0195     return;
0196   }
0197 
0198   if (!mask)
0199   {
0200     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_z - invalid mask histogram" << std::endl;
0201     return;
0202   }
0203 
0204   // loop over phi bins
0205   for (int ir = 0; ir < source->GetNbinsY(); ++ir)
0206   {
0207     for (int ip = 0; ip < source->GetNbinsX(); ++ip)
0208     {
0209       // handle pad plane on the positive side (zmax)
0210       if (side & Side_positive)
0211       {
0212         // find first non empty bin, starting from last zaxis bin
0213         int iz_min = -1;
0214         for (int iz = source->GetNbinsZ() - 1; iz >= 0; --iz)
0215         {
0216           const bool in_range = mask->GetBinContent(ip + 1, ir + 1, iz + 1) > 0;
0217           if (in_range)
0218           {
0219             iz_min = iz;
0220             break;
0221           }
0222         }
0223 
0224         // check bin validity
0225         if (iz_min < 0)
0226         {
0227           continue;
0228         }
0229 
0230         // loop over empty bins and do the interpolation
0231         for (int iz = iz_min + 1; iz < source->GetNbinsZ(); ++iz)
0232         {
0233           const double z_min = source->GetZaxis()->GetBinCenter(iz_min + 1);
0234           const double z_max = source->GetZaxis()->GetXmax();
0235           const double z = source->GetZaxis()->GetBinCenter(iz + 1);
0236           const double alpha_min = (z_max - z) / (z_max - z_min);
0237           const double alpha_max = (z - z_min) / (z_max - z_min);
0238 
0239           const double distortion_min = source->GetBinContent(ip + 1, ir + 1, iz_min + 1);
0240           const double distortion_max = 0;
0241           const double distortion = alpha_min * distortion_min + alpha_max * distortion_max;
0242 
0243           const double error_min = source->GetBinError(ip + 1, ir + 1, iz_min + 1);
0244           const double error_max = 0;
0245           const double error = std::sqrt(square(alpha_min * error_min) + square(alpha_max * error_max));
0246 
0247           source->SetBinContent(ip + 1, ir + 1, iz + 1, distortion);
0248           source->SetBinError(ip + 1, ir + 1, iz + 1, error);
0249         }
0250       }
0251 
0252       // handle pad plane on the negative side (zmin
0253       if (side & Side_negative)
0254       {
0255         // find first non empty bin, starting from first zaxis bin
0256         int iz_max = -1;
0257         for (int iz = 0; iz < source->GetNbinsZ(); ++iz)
0258         {
0259           const bool in_range = mask->GetBinContent(ip + 1, ir + 1, iz + 1) > 0;
0260           if (in_range)
0261           {
0262             iz_max = iz;
0263             break;
0264           }
0265         }
0266 
0267         // check bin validity
0268         if (iz_max < 0)
0269         {
0270           continue;
0271         }
0272 
0273         // loop over empty bins and do the interpolation
0274         for (int iz = 0; iz < iz_max; ++iz)
0275         {
0276           const double z_min = source->GetZaxis()->GetXmin();
0277           const double z_max = source->GetZaxis()->GetBinCenter(iz_max + 1);
0278           const double z = source->GetZaxis()->GetBinCenter(iz + 1);
0279           const double alpha_min = (z_max - z) / (z_max - z_min);
0280           const double alpha_max = (z - z_min) / (z_max - z_min);
0281 
0282           const double distortion_min = 0;
0283           const double distortion_max = source->GetBinContent(ip + 1, ir + 1, iz_max + 1);
0284           const double distortion = alpha_min * distortion_min + alpha_max * distortion_max;
0285 
0286           const double error_min = 0;
0287           const double error_max = source->GetBinError(ip + 1, ir + 1, iz_max + 1);
0288           const double error = std::sqrt(square(alpha_min * error_min) + square(alpha_max * error_max));
0289 
0290           source->SetBinContent(ip + 1, ir + 1, iz + 1, distortion);
0291           source->SetBinError(ip + 1, ir + 1, iz + 1, error);
0292         }
0293       }
0294     }
0295   }
0296 }
0297 
0298 //____________________________________________________________________________________
0299 void TpcSpaceChargeReconstructionHelper::extrapolate_phi1(TH3* source, const TH2* source_cm, const TH3* mask)
0300 {
0301   if (!source)
0302   {
0303     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_phi1 - invalid source histogram" << std::endl;
0304     return;
0305   }
0306 
0307   if (!mask)
0308   {
0309     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_phi1 - invalid mask histogram" << std::endl;
0310     return;
0311   }
0312 
0313   // loop over phi bins
0314   for (int ip = 0; ip < source->GetNbinsX(); ++ip)
0315   {
0316 
0317     // get phi
0318     const double phi = source->GetXaxis()->GetBinCenter(ip + 1);
0319     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_phi1 -"
0320       << " source: " << source->GetName()
0321       << " ip: " << ip
0322       << " phi: " << phi
0323       << std::endl;
0324 
0325     for (int ir = 0; ir < source->GetNbinsY(); ++ir)
0326     {
0327 
0328       // get radius
0329       const double r = source->GetYaxis()->GetBinCenter(ir + 1);
0330       std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_phi1 -"
0331         << " source: " << source->GetName()
0332         << " ir: " << ir
0333         << " r: " << r
0334         << std::endl;
0335 
0336       for (int iz = 0; iz < source->GetNbinsZ(); ++iz)
0337       {
0338         // do nothing if in TPOT acceptance
0339         bool in_range = mask->GetBinContent(ip + 1, ir + 1, iz + 1) > 0;
0340         if (in_range)
0341         {
0342           continue;
0343         }
0344 
0345         // phi not in TPOT range, rotate by steps of 2pi/12 (= one TPC sector) until found in range
0346         static constexpr int n_sectors = 12;
0347         for (int sector = 1; sector < n_sectors; ++sector)
0348         {
0349           // get ref phi
0350           double phi_ref = phi + 2. * M_PI * sector / n_sectors;
0351           while (phi_ref >= 2 * M_PI)
0352           {
0353             phi_ref -= 2 * M_PI;
0354           };
0355 
0356           int ip_ref = source->GetXaxis()->FindBin(phi_ref) - 1;
0357           in_range = mask->GetBinContent(ip_ref + 1, ir + 1, iz + 1) > 0;
0358           if (!in_range)
0359           {
0360             continue;
0361           }
0362 
0363           // get normalization factor from CM histograms
0364           double scale = 1;
0365           if (source_cm)
0366           {
0367 
0368             auto ip_cm = source_cm->GetXaxis()->FindBin(phi);
0369             auto ip_ref_cm = source_cm->GetXaxis()->FindBin(phi_ref);
0370             auto ir_cm = source_cm->GetYaxis()->FindBin( r );
0371 
0372             // check out of bound
0373             if( ( ip_cm == 0 || ip_cm > source_cm->GetNbinsX() ) ||
0374               ( ir==0 || ir > source_cm->GetNbinsY() ) )
0375             {
0376 
0377               scale = 1;
0378 
0379             } else if( ir_cm==1 || ir_cm == source_cm->GetNbinsY() ) {
0380 
0381               // if first or last bin is used, interpolate will break. Need to use the bin content
0382               const double distortion_local = source_cm->GetBinContent(ip_cm, ir_cm);
0383               const double distortion_ref = source_cm->GetBinContent(ip_ref_cm, ir_cm);
0384               scale = distortion_local / distortion_ref;
0385 
0386             } else if( ip_cm==1 || ip_cm == source_cm->GetNbinsX() ) {
0387 
0388               // if first or last bin is used, interpolate will break. Need to use the bin content
0389               const double distortion_local = source_cm->GetBinContent(ip_cm, ir_cm);
0390               const double distortion_ref = source_cm->GetBinContent(ip_ref_cm, ir_cm);
0391               scale = distortion_local / distortion_ref;
0392 
0393             } else {
0394 
0395               // use interpolation
0396               const double distortion_local = source_cm->Interpolate(phi, r);
0397               const double distortion_ref = source_cm->Interpolate(phi_ref, r);
0398               scale = distortion_local / distortion_ref;
0399 
0400             }
0401 
0402             // check imposible values
0403             if( std::isnan(scale) ) { scale = 1; }
0404           }
0405 
0406           // update content and error
0407           const double distortion = scale * source->GetBinContent(ip_ref + 1, ir + 1, iz + 1);
0408           const double error = scale * source->GetBinError(ip_ref + 1, ir + 1, iz + 1);
0409           source->SetBinContent(ip + 1, ir + 1, iz + 1, distortion);
0410           source->SetBinError(ip + 1, ir + 1, iz + 1, error);
0411         }
0412       }
0413     }
0414   }
0415 }
0416 
0417 //_______________________________________________
0418 void TpcSpaceChargeReconstructionHelper::extrapolate_phi2(TH3* source, const TH3* mask)
0419 {
0420   if (!source)
0421   {
0422     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_phi2 - invalid source histogram" << std::endl;
0423     return;
0424   }
0425 
0426   if (!mask)
0427   {
0428     std::cout << "TpcSpaceChargeReconstructionHelper::extrapolate_phi2 - invalid mask histogram" << std::endl;
0429     return;
0430   }
0431 
0432   // loop over phi bins
0433   for (int ir = 0; ir < source->GetNbinsY(); ++ir)
0434   {
0435     for (int iz = 0; iz < source->GetNbinsZ(); ++iz)
0436     {
0437       int ip_min = -1;
0438       for (int ip = 0; ip < source->GetNbinsX(); ++ip)
0439       {
0440         bool in_range = mask->GetBinContent(ip + 1, ir + 1, iz + 1) > 0;
0441         if (in_range)
0442         {
0443           ip_min = ip;
0444           continue;
0445         }
0446 
0447         // ip is not in range. Check if ip_min was set. If not, skip this bin. If yes, find next bin in range
0448         if (ip_min < 0)
0449         {
0450           continue;
0451         }
0452 
0453         int ip_max = ip + 1;
0454         for (; ip_max < source->GetNbinsX(); ++ip_max)
0455         {
0456           in_range = mask->GetBinContent(ip_max + 1, ir + 1, iz + 1) > 0;
0457           if (in_range)
0458           {
0459             break;
0460           }
0461         }
0462 
0463         // check that a valid bin was found
0464         if (ip_max == source->GetNbinsX())
0465         {
0466           continue;
0467         }
0468 
0469         // do the interpolation
0470         const double phi_min = source->GetXaxis()->GetBinCenter(ip_min + 1);
0471         const double phi_max = source->GetXaxis()->GetBinCenter(ip_max + 1);
0472         const double phi = source->GetXaxis()->GetBinCenter(ip + 1);
0473 
0474         const double alpha_min = (phi_max - phi) / (phi_max - phi_min);
0475         const double alpha_max = (phi - phi_min) / (phi_max - phi_min);
0476 
0477         const double distortion_min = source->GetBinContent(ip_min + 1, ir + 1, iz + 1);
0478         const double distortion_max = source->GetBinContent(ip_max + 1, ir + 1, iz + 1);
0479         const double distortion = alpha_min * distortion_min + alpha_max * distortion_max;
0480 
0481         const double error_min = source->GetBinError(ip_min + 1, ir + 1, iz + 1);
0482         const double error_max = source->GetBinError(ip_max + 1, ir + 1, iz + 1);
0483         const double error = std::sqrt(square(alpha_min * error_min) + square(alpha_max * error_max));
0484 
0485         source->SetBinContent(ip + 1, ir + 1, iz + 1, distortion);
0486         source->SetBinError(ip + 1, ir + 1, iz + 1, error);
0487       }
0488     }
0489   }
0490 }
0491 
0492 //_______________________________________________
0493 std::tuple<TH3*, TH3*> TpcSpaceChargeReconstructionHelper::split(const TH3* source)
0494 {
0495   if (!source)
0496   {
0497     return std::make_tuple<TH3*, TH3*>(nullptr, nullptr);
0498   }
0499 
0500   auto xaxis = source->GetXaxis();
0501   auto yaxis = source->GetYaxis();
0502   auto zaxis = source->GetZaxis();
0503   auto ibin = zaxis->FindBin((double) 0);
0504 
0505   // create histograms
0506   const TString name(source->GetName());
0507   auto hneg = new TH3F(
0508       name + "_negz", name + "_negz",
0509       xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax(),
0510       yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax(),
0511       ibin - 1, zaxis->GetXmin(), zaxis->GetBinUpEdge(ibin - 1));
0512 
0513   auto hpos = new TH3F(
0514       name + "_posz", name + "_posz",
0515       xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax(),
0516       yaxis->GetNbins(), yaxis->GetXmin(), yaxis->GetXmax(),
0517       zaxis->GetNbins() - (ibin - 1), zaxis->GetBinLowEdge(ibin), zaxis->GetXmax());
0518 
0519   // copy content and errors
0520   for (int ix = 0; ix < xaxis->GetNbins(); ++ix)
0521   {
0522     for (int iy = 0; iy < yaxis->GetNbins(); ++iy)
0523     {
0524       for (int iz = 0; iz < zaxis->GetNbins(); ++iz)
0525       {
0526         const auto content = source->GetBinContent(ix + 1, iy + 1, iz + 1);
0527         const auto error = source->GetBinError(ix + 1, iy + 1, iz + 1);
0528 
0529         if (iz < ibin - 1)
0530         {
0531           hneg->SetBinContent(ix + 1, iy + 1, iz + 1, content);
0532           hneg->SetBinError(ix + 1, iy + 1, iz + 1, error);
0533         }
0534         else
0535         {
0536           hpos->SetBinContent(ix + 1, iy + 1, iz - (ibin - 1) + 1, content);
0537           hpos->SetBinError(ix + 1, iy + 1, iz - (ibin - 1) + 1, error);
0538         }
0539       }
0540     }
0541   }
0542 
0543   // also copy axis titles
0544   for (const auto h : {hneg, hpos})
0545   {
0546     h->GetXaxis()->SetTitle(source->GetXaxis()->GetTitle());
0547     h->GetYaxis()->SetTitle(source->GetYaxis()->GetTitle());
0548     h->GetZaxis()->SetTitle(source->GetZaxis()->GetTitle());
0549   }
0550 
0551   return std::make_tuple(hneg, hpos);
0552 }
0553 
0554 //___________________________________________________________________________
0555 TH3* TpcSpaceChargeReconstructionHelper::add_guarding_bins(const TH3* source, const TString& name)
0556 {
0557   std::array<int, 3> bins{};
0558   std::array<double, 3> x_min{};
0559   std::array<double, 3> x_max{};
0560 
0561   int index = 0;
0562   for (const auto axis : {source->GetXaxis(), source->GetYaxis(), source->GetZaxis()})
0563   {
0564     // calculate bin width
0565     const auto bin_width = (axis->GetXmax() - axis->GetXmin()) / axis->GetNbins();
0566 
0567     // increase the number of bins by two
0568     bins[index] = axis->GetNbins() + 2;
0569 
0570     // update axis limits accordingly
0571     x_min[index] = axis->GetXmin() - bin_width;
0572     x_max[index] = axis->GetXmax() + bin_width;
0573     ++index;
0574   }
0575 
0576   // create new histogram
0577   auto hout = new TH3F(name, name,
0578                        bins[0], x_min[0], x_max[0],
0579                        bins[1], x_min[1], x_max[1],
0580                        bins[2], x_min[2], x_max[2]);
0581 
0582   // update axis legend
0583   hout->GetXaxis()->SetTitle(source->GetXaxis()->GetTitle());
0584   hout->GetYaxis()->SetTitle(source->GetYaxis()->GetTitle());
0585   hout->GetZaxis()->SetTitle(source->GetZaxis()->GetTitle());
0586 
0587   // copy content
0588   const auto phibins = source->GetXaxis()->GetNbins();
0589   const auto rbins = source->GetYaxis()->GetNbins();
0590   const auto zbins = source->GetZaxis()->GetNbins();
0591 
0592   // fill center
0593   for (int iphi = 0; iphi < phibins; ++iphi)
0594   {
0595     for (int ir = 0; ir < rbins; ++ir)
0596     {
0597       for (int iz = 0; iz < zbins; ++iz)
0598       {
0599         hout->SetBinContent(iphi + 2, ir + 2, iz + 2, source->GetBinContent(iphi + 1, ir + 1, iz + 1));
0600         hout->SetBinError(iphi + 2, ir + 2, iz + 2, source->GetBinError(iphi + 1, ir + 1, iz + 1));
0601       }
0602     }
0603   }
0604 
0605   // fill guarding phi bins
0606   /*
0607    * we use 2pi periodicity to do that:
0608    * - last valid bin is copied to first guarding bin;
0609    * - first valid bin is copied to last guarding bin
0610    */
0611   for (int ir = 0; ir < rbins + 2; ++ir)
0612   {
0613     for (int iz = 0; iz < zbins + 2; ++iz)
0614     {
0615       // copy last bin to first guarding bin
0616       hout->SetBinContent(1, ir + 1, iz + 1, hout->GetBinContent(phibins + 1, ir + 1, iz + 1));
0617       hout->SetBinError(1, ir + 1, iz + 1, hout->GetBinError(phibins + 1, ir + 1, iz + 1));
0618 
0619       // copy first bin to last guarding bin
0620       hout->SetBinContent(phibins + 2, ir + 1, iz + 1, hout->GetBinContent(2, ir + 1, iz + 1));
0621       hout->SetBinError(phibins + 2, ir + 1, iz + 1, hout->GetBinError(2, ir + 1, iz + 1));
0622     }
0623   }
0624 
0625   // fill guarding r bins
0626   for (int iphi = 0; iphi < phibins + 2; ++iphi)
0627   {
0628     for (int iz = 0; iz < zbins + 2; ++iz)
0629     {
0630       hout->SetBinContent(iphi + 1, 1, iz + 1, hout->GetBinContent(iphi + 1, 2, iz + 1));
0631       hout->SetBinError(iphi + 1, 1, iz + 1, hout->GetBinError(iphi + 1, 2, iz + 1));
0632 
0633       hout->SetBinContent(iphi + 1, rbins + 2, iz + 1, hout->GetBinContent(iphi + 1, rbins + 1, iz + 1));
0634       hout->SetBinError(iphi + 1, rbins + 2, iz + 1, hout->GetBinError(iphi + 1, rbins + 1, iz + 1));
0635     }
0636   }
0637 
0638   // fill guarding z bins
0639   for (int iphi = 0; iphi < phibins + 2; ++iphi)
0640   {
0641     for (int ir = 0; ir < rbins + 2; ++ir)
0642     {
0643       hout->SetBinContent(iphi + 1, ir + 1, 1, hout->GetBinContent(iphi + 1, ir + 1, 2));
0644       hout->SetBinError(iphi + 1, ir + 1, 1, hout->GetBinError(iphi + 1, ir + 1, 2));
0645 
0646       hout->SetBinContent(iphi + 1, ir + 1, zbins + 2, hout->GetBinContent(iphi + 1, ir + 1, zbins + 1));
0647       hout->SetBinError(iphi + 1, ir + 1, zbins + 2, hout->GetBinError(iphi + 1, ir + 1, zbins + 1));
0648     }
0649   }
0650 
0651   return hout;
0652 }
0653 
0654 //___________________________________________________________________________________________________
0655 void TpcSpaceChargeReconstructionHelper::set_phi_range_central( const range_t& range)
0656 { phi_range_central = transform_range(range); }
0657 
0658 //___________________________________________________________________________________________________
0659 void TpcSpaceChargeReconstructionHelper::set_phi_range_east( const range_t& range)
0660 { phi_range_east = transform_range(range); }
0661 
0662 //___________________________________________________________________________________________________
0663 void TpcSpaceChargeReconstructionHelper::set_phi_range_west( const range_t& range)
0664 { phi_range_west = transform_range(range); }