packages/PHGarfield/PHGarfieldRossegger.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef PHGARFIELD_PHGARFIELDROSSEGGER_H
0004 #define PHGARFIELD_PHGARFIELDROSSEGGER_H
0005
0006 #include <fun4all/SubsysReco.h>
0007
0008 #include <array>
0009 #include <string>
0010 #include <utility>
0011 #include <vector>
0012
0013 class PHCompositeNode;
0014
0015 class PHGarfieldRossegger : public SubsysReco
0016 {
0017  public:
0018   explicit PHGarfieldRossegger(const std::string& name = "PHGarfieldRossegger");
0019   ~PHGarfieldRossegger() override = default;
0020
0021   int Init(PHCompositeNode*) override;
0022
0023   int InitRun(PHCompositeNode*) override;
0024   int process_event(PHCompositeNode*) override;
0025
0026   void setGeometryCm(double inner_radius, double outer_radius, double half_length);
0027   void setSourceRadiusCm(double min_radius, double max_radius);
0028   void setDensity(double reference_density_nC_per_m3, double k_eff, double radial_power_alpha);
0029   void setPhiModulation(double m1_amplitude, double m1_phase, double m12_amplitude, double m12_phase);
0030   void setSourceGrid(unsigned int nr, unsigned int nphi, unsigned int nz);
0031   void setObservationGrid(unsigned int nr, unsigned int nphi, unsigned int nz);
0032   void setModeTruncation(unsigned int m_phi_max, unsigned int n_radial_modes, unsigned int n_longitudinal_modes);
0033   void setAutoAxisymmetric(bool value) { m_autoAxisymmetric = value; }
0034   void setTpcSide(unsigned int side) { m_tpcSide = side; }
0035
0036   void setUseRealTpcSourceGeometry(bool value) { m_useRealTpcSourceGeometry = value; }
0037   void setPadPlacementFile(const std::string& filename) { m_padPlacementFile = filename; }
0038   void setGainMapFile(const std::string& filename) { m_gainMapFile = filename; }
0039   void setGainHistograms(const std::string& side0_histogram, const std::string& side1_histogram);
0040   void setDivideChargeByGain(bool value) { m_divideChargeByGain = value; }
0041   void setNormalizeGainWeightedTotal(bool value) { m_normalizeGainWeightedTotal = value; }
0042
0043   void setUseDensityMap(bool value) { m_useDensityMap = value; }
0044   void setDensityMapFile(const std::string& filename, const std::string& side0_histogram, const std::string& side1_histogram);
0045   void setNormalizeDensityMap(bool value) { m_normalizeDensityMap = value; }
0046
0047   void setUseFrameChargeModel(bool value) { m_useFrameChargeModel = value; }
0048   void setFrameReferencePhi(double value) { m_frameReferencePhi = value; }
0049   void setFrameGeometryFile(const std::string& filename) { m_frameGeometryFile = filename; m_framePolygonsLoaded = false; }
0050   enum class FrameChargeWeighting
0051   {
0052     EqualChargePerPiece,
0053     ProportionalToArea
0054   };
0055   void setFrameBoundaryPotential(double value) { m_frameBoundaryPotential = value; }
0056   void setFrameChargeWeighting(FrameChargeWeighting mode) { m_frameChargeWeighting = mode; }
0057   void setFrameEzScale(double value) { m_frameEzScale = value; }
0058
0059   // Split observation radial bins across condor jobs. Histograms keep full
0060   // binning, so PART files can be merged with hadd.
0061   void setRadialJob(unsigned int job_index, unsigned int n_jobs);
0062
0063   void setGarfieldOutputFile(const std::string& filename) { m_garfieldOutputFile = filename; }
0064   void setField3DOutputFile(const std::string& filename) { m_field3DOutputFile = filename; }
0065   void setPHGarfieldField3DOutputFiles(const std::string& side0_filename, const std::string& side1_filename);
0066   void setWriteField3D(bool value) { m_writeField3D = value; }
0067   void setWritePHGarfieldField3D(bool value) { m_writePHGarfieldField3D = value; }
0068   void setVerifyOutput(bool value) { m_verifyOutput = value; }
0069
0070  private:
0071   struct SourceGrid
0072   {
0073     std::vector<double> r_edges_m;
0074     std::vector<double> r_centers_m;
0075     std::vector<int> module_index;
0076     std::vector<int> layer_index;
0077     std::vector<bool> is_antenna;
0078   };
0079
0080   struct PadInterval
0081   {
0082     double start_edge{0.0};
0083     double angular_width{0.0};
0084     int direction{1};
0085   };
0086
0087   using PadGeometry = std::array<std::array<std::array<PadInterval, 12>, 2>, 3>;
0088   using GainMap = std::array<std::array<std::array<double, 48>, 12>, 2>;
0089
0090   struct Point2D
0091   {
0092     double x_cm{0.0};
0093     double y_cm{0.0};
0094   };
0095
0096   struct FramePolygon
0097   {
0098     std::vector<Point2D> inner;
0099     std::vector<Point2D> outer;
0100   };
0101
0102   using FramePolygons = std::array<FramePolygon, 3>;
0103
0104   struct FrameBoundaryPattern
0105   {
0106     std::vector<double> geometry_fraction;
0107     std::vector<double> weight;
0108     std::vector<double> boundary_potential;
0109   };
0110
0111   int calculate();
0112   bool validateConfig() const;
0113   unsigned int effectivePhiMax() const;
0114   std::pair<unsigned int, unsigned int> radialRange(unsigned int nr_obs) const;
0115
0116   double radialBoundaryFunction(double kval, unsigned int mode_m) const;
0117   double radialBasis(double kval, unsigned int mode_m, double radius_m) const;
0118   double radialBasisDerivative(double kval, unsigned int mode_m, double radius_m) const;
0119   std::vector<double> findRadialRoots(unsigned int mode_m, unsigned int n_roots) const;
0120   std::vector<double> legendreRootsAndWeights(unsigned int n_points, std::vector<double>& weights) const;
0121   SourceGrid buildSourceGrid() const;
0122   PadGeometry parsePadPlacementFile() const;
0123   GainMap readGainMap() const;
0124   std::vector<double> makeChargeDensity(const SourceGrid& source_grid,
0125                                         const std::vector<double>& phi_source_edges,
0126                                         const std::vector<double>& z_source_edges_m,
0127                                         unsigned int side) const;
0128   void loadFramePolygons() const;
0129   static bool pointInPolygon(const std::vector<Point2D>& polygon, double x_cm, double y_cm);
0130   static double polygonAreaCm2(const std::vector<Point2D>& polygon);
0131   bool pointInFrameGeometry(const FramePolygon& frame, double r_cm, double phi_rel) const;
0132   double frameAreaM2(const FramePolygon& frame) const;
0133   FrameBoundaryPattern makeFrameBoundaryPattern(const SourceGrid& source_grid,
0134                                                 const std::vector<double>& phi_source_edges,
0135                                                 unsigned int side) const;
0136
0137   void writeGarfieldRootFile(const std::vector<double>& r_edges_m,
0138                              const std::vector<double>& phi_edges,
0139                              const std::vector<double>& z_edges_m,
0140                              const std::vector<double>& er,
0141                              const std::vector<double>& ephi,
0142                              const std::vector<double>& ez,
0143                              unsigned int r_begin,
0144                              unsigned int r_end) const;
0145
0146   void writeField3DRootFile(const std::vector<double>& r_source_edges_m,
0147                             const std::vector<double>& r_obs_edges_m,
0148                             const std::vector<double>& phi_source_edges,
0149                             const std::vector<double>& phi_obs_edges,
0150                             const std::vector<double>& z_source_edges_m,
0151                             const std::vector<double>& z_obs_edges_m,
0152                             const std::vector<double>& rho,
0153                             const std::vector<double>& potential,
0154                             const std::vector<double>& er,
0155                             const std::vector<double>& ephi,
0156                             const std::vector<double>& ez,
0157                             unsigned int r_begin,
0158                             unsigned int r_end) const;
0159
0160   void writeFrameBoundaryField3DRootFile(const std::vector<double>& r_source_edges_m,
0161                                          const std::vector<double>& r_obs_edges_m,
0162                                          const std::vector<double>& phi_source_edges,
0163                                          const std::vector<double>& phi_obs_edges,
0164                                          const std::vector<double>& z_obs_edges_m,
0165                                          const FrameBoundaryPattern& frame_pattern,
0166                                          const std::vector<double>& potential,
0167                                          const std::vector<double>& er,
0168                                          const std::vector<double>& ephi,
0169                                          const std::vector<double>& ez,
0170                                          unsigned int r_begin,
0171                                          unsigned int r_end) const;
0172
0173   void writePHGarfieldField3DRootFile(const std::vector<double>& r_obs_edges_m,
0174                                       const std::vector<double>& phi_obs_edges,
0175                                       const std::vector<double>& z_obs_edges_m,
0176                                       const std::vector<double>& potential,
0177                       const std::vector<double>& er,
0178                                       const std::vector<double>& ephi,
0179                                       const std::vector<double>& ez,
0180                                       unsigned int side,
0181                                       unsigned int r_begin,
0182                                       unsigned int r_end) const;
0183
0184   bool verifyOutput() const;
0185
0186   double m_aCm{21.6};
0187   double m_bCm{76.4};
0188   double m_lCm{102.325};
0189   double m_sourceRMinCm{22.8};
0190   double m_sourceRMaxCm{75.43};
0191
0192   double m_rhoReferenceNCPerM3{20.0};
0193   double m_kEff{1.0};
0194   double m_radialPowerAlpha{1.4};
0195   double m_m1Amplitude{0.0};
0196   double m_m1Phase{0.0};
0197   double m_m12Amplitude{0.0};
0198   double m_m12Phase{0.0};
0199
0200   unsigned int m_tpcSide{0};
0201   unsigned int m_nrSource{18};
0202   unsigned int m_nphiSource{24};
0203   unsigned int m_nzSource{20};
0204   unsigned int m_nrObs{24};
0205   unsigned int m_nphiObs{24};
0206   unsigned int m_nzObs{28};
0207   unsigned int m_mPhiMax{4};
0208   unsigned int m_nRadialModes{18};
0209   unsigned int m_nLongitudinalModes{22};
0210   bool m_autoAxisymmetric{true};
0211   bool m_useRealTpcSourceGeometry{false};
0212   bool m_useFrameChargeModel{false};
0213   double m_frameReferencePhi{0.0};
0214   std::string m_frameGeometryFile{"tpc_frame_geometry.csv"};
0215   mutable FramePolygons m_framePolygons;
0216   mutable bool m_framePolygonsLoaded{false};
0217   double m_frameBoundaryPotential{1.0};
0218   FrameChargeWeighting m_frameChargeWeighting{FrameChargeWeighting::ProportionalToArea};
0219   double m_frameEzScale{1.0};
0220   std::string m_padPlacementFile{"TPC_pad_placement.txt"};
0221   std::string m_gainMapFile{"layer_gain_79513_Mariia_side01.root"};
0222   std::array<std::string, 2> m_gainHistograms{{"hGainMap_side0_South", "hGainMap_side1_North"}};
0223   bool m_divideChargeByGain{true};
0224   bool m_normalizeGainWeightedTotal{true};
0225
0226   bool m_useDensityMap{false};
0227   bool m_normalizeDensityMap{true};
0228   std::string m_densityMapFile;
0229   std::array<std::string, 2> m_densityMapHistograms{{"h_ibf_final_rphi_side0", "h_ibf_final_rphi_side1"}};
0230
0231   unsigned int m_jobIndex{0};
0232   unsigned int m_nJobs{1};
0233
0234   std::string m_garfieldOutputFile{"sphenix_rossegger_garfield_field_1p4.root"};
0235   std::string m_field3DOutputFile{"sphenix_rossegger_fields_3d.root"};
0236   std::array<std::string, 2> m_phGarfieldField3DOutputFiles{{"sphenix_3d_ibf_field_side0_South_v1.root", "sphenix_3d_ibf_field_side1_North_v1.root"}};
0237   bool m_writeField3D{false};
0238   bool m_writePHGarfieldField3D{false};
0239   bool m_verifyOutput{true};
0240   bool m_done{false};
0241 };
0242
0243 #endif