g4simulation/g4tpc/PHG4TpcCentralMembrane.h

0001 #ifndef PHG4TPCCENTRALMEMBRANE_H
0002 #define PHG4TPCCENTRALMEMBRANE_H
0003
0004 #include <phparameter/PHParameterInterface.h>
0005
0006 #include <fun4all/SubsysReco.h>
0007
0008 #include <array>
0009 #include <cmath>
0010 #include <limits>
0011 #include <string>  // for string
0012 #include <vector>
0013
0014 class PHCompositeNode;
0015 class PHG4Hit;
0016
0017 // all distances in mm, all angles in rad
0018 // class that generates stripes and dummy hit coordinates
0019 // stripes have width of one mm, length of one pad width, and are centered in middle of sector gaps
0020
0021 class PHG4TpcCentralMembrane : public SubsysReco, public PHParameterInterface
0022 {
0023  public:
0024   /// constructor
0025   PHG4TpcCentralMembrane(const std::string& name = "PHG4TpcCentralMembrane");
0026
0027   /// destructor
0028   ~PHG4TpcCentralMembrane() override;
0029
0030   /// run initialization
0031   int InitRun(PHCompositeNode*) override;
0032
0033   /// per event processing
0034   int process_event(PHCompositeNode*) override;
0035
0036   /// default parameters
0037   void SetDefaultParameters() override;
0038
0039   /// detector name
0040   void Detector(const std::string& d)
0041   {
0042     detector = d;
0043     hitnodename = "G4HIT_" + d;
0044   }
0045
0046   /// check if coords are in a stripe
0047   int getSearchResult(double xcheck, double ycheck) const;
0048
0049   int getStripeID(double xcheck, double ycheck) const;
0050
0051   /// adjust central membrane hits delay with respect to trigger time
0052   void setCentralMembraneDelay(int ns) { m_centralMembraneDelay = ns; };
0053
0054   /// set modulo for events in which to generate CM hits
0055   void setCentralMembraneEventModulo(int mod) { m_eventModulo = mod; };
0056
0057  private:
0058   /// detector name
0059   std::string detector = "TPC";
0060
0061   /// g4hitnode name
0062   std::string hitnodename = "G4HIT_TPC";
0063   std::vector<PHG4Hit*> PHG4Hits;
0064
0065   int m_eventModulo = 10;
0066   int m_eventNum = 0;
0067
0068   static constexpr double mm = 1.0;
0069   static constexpr double cm = 10.0;
0070
0071   /// inner radius of CM
0072   static constexpr double begin_CM = 221.4019814 * mm;
0073
0074   /// outer radius of CM
0075   static constexpr double end_CM = 759.2138 * mm;
0076
0077   static constexpr int nRadii = 8;
0078   static constexpr int nStripes_R1 = 6;
0079   static constexpr int nStripes_R2 = 8;
0080   static constexpr int nStripes_R3 = 12;
0081
0082   static constexpr int nPads_R1 = 6 * 16;
0083   static constexpr int nPads_R2 = 8 * 16;
0084   static constexpr int nPads_R3 = 12 * 16;
0085
0086   //
0087   static constexpr double padfrac_R1 = 0.5 * 5.59106385 * mm;
0088   static constexpr double padfrac_R2 = 0.5 * 10.13836283 * mm;
0089   static constexpr double padfrac_R3 = 0.5 * 10.90189537 * mm;
0090
0091   /// radius of arc on end of a stripe
0092   static constexpr double arc_r = 0.5 * mm;
0093
0094   /// stripe radii
0095   static constexpr std::array<double, nRadii> R1_e = {{227.0902789 * mm, 238.4100043 * mm, 249.7297296 * mm, 261.049455 * mm, 272.3691804 * mm, 283.6889058 * mm, 295.0086312 * mm, 306.3283566 * mm}};
0096   static constexpr std::array<double, nRadii> R1 = {{317.648082 * mm, 328.9678074 * mm, 340.2875328 * mm, 351.6072582 * mm, 362.9269836 * mm, 374.246709 * mm, 385.5664344 * mm, 396.8861597 * mm}};
0097   static constexpr std::array<double, nRadii> R2 = {{421.705532 * mm, 442.119258 * mm, 462.532984 * mm, 482.9467608 * mm, 503.36069 * mm, 523.774416 * mm, 544.188015 * mm, 564.601868 * mm}};
0098   static constexpr std::array<double, nRadii> R3 = {{594.6048725 * mm, 616.545823 * mm, 638.4867738 * mm, 660.4277246 * mm, 682.3686754 * mm, 704.3096262 * mm, 726.250577 * mm, 748.1915277 * mm}};
0099
0100   double str_width_R1_e[nStripes_R1][nRadii] = {};
0101   double str_width_R1[nStripes_R1][nRadii] = {};
0102   double str_width_R2[nStripes_R2][nRadii] = {};
0103   double str_width_R3[nStripes_R3][nRadii] = {};
0104
0105   static constexpr std::array<double, nRadii> widthmod_R1_e = {{1.493, 1.398, 1.334, 1.284, 1.243, 1.208, 1.178, 1.152}};
0106   static constexpr std::array<double, nRadii> widthmod_R1 = {{1.129, 1.109, 1.091, 1.076, 1.062, 1.050, 1.040, 1.030}};
0107   static constexpr std::array<double, nRadii> widthmod_R2 = {{1.015, 1.007, 1.002, 1.000, 1.001, 1.006, 1.013, 1.023}};
0108   static constexpr std::array<double, nRadii> widthmod_R3 = {{1.044, 1.064, 1.087, 1.115, 1.147, 1.186, 1.232, 1.288}};
0109
0110   std::array<double, nRadii> spacing_R1_e{};
0111   std::array<double, nRadii> spacing_R1{};
0112   std::array<double, nRadii> spacing_R2{};
0113   std::array<double, nRadii> spacing_R3{};
0114
0115   // bottom left - 1a
0116   double x1a_R1_e[nStripes_R1][nRadii]{};
0117   double y1a_R1_e[nStripes_R1][nRadii]{};
0118   double x1a_R1[nStripes_R1][nRadii]{};
0119   double y1a_R1[nStripes_R1][nRadii]{};
0120   double x1a_R2[nStripes_R2][nRadii]{};
0121   double y1a_R2[nStripes_R2][nRadii]{};
0122   double x1a_R3[nStripes_R3][nRadii]{};
0123   double y1a_R3[nStripes_R3][nRadii]{};
0124
0125   // bottom right - 1b
0126   double x1b_R1_e[nStripes_R1][nRadii]{};
0127   double y1b_R1_e[nStripes_R1][nRadii]{};
0128   double x1b_R1[nStripes_R1][nRadii]{};
0129   double y1b_R1[nStripes_R1][nRadii]{};
0130   double x1b_R2[nStripes_R2][nRadii]{};
0131   double y1b_R2[nStripes_R2][nRadii]{};
0132   double x1b_R3[nStripes_R3][nRadii]{};
0133   double y1b_R3[nStripes_R3][nRadii]{};
0134
0135   // top left - 2a
0136   double x2a_R1_e[nStripes_R1][nRadii]{};
0137   double y2a_R1_e[nStripes_R1][nRadii]{};
0138   double x2a_R1[nStripes_R1][nRadii]{};
0139   double y2a_R1[nStripes_R1][nRadii]{};
0140   double x2a_R2[nStripes_R2][nRadii]{};
0141   double y2a_R2[nStripes_R2][nRadii]{};
0142   double x2a_R3[nStripes_R3][nRadii]{};
0143   double y2a_R3[nStripes_R3][nRadii]{};
0144
0145   // top right - 2b
0146   double x2b_R1_e[nStripes_R1][nRadii]{};
0147   double y2b_R1_e[nStripes_R1][nRadii]{};
0148   double x2b_R1[nStripes_R1][nRadii]{};
0149   double y2b_R1[nStripes_R1][nRadii]{};
0150   double x2b_R2[nStripes_R2][nRadii]{};
0151   double y2b_R2[nStripes_R2][nRadii]{};
0152   double x2b_R3[nStripes_R3][nRadii]{};
0153   double y2b_R3[nStripes_R3][nRadii]{};
0154
0155   // left midpoint - 3a
0156   double x3a_R1_e[nStripes_R1][nRadii]{};
0157   double y3a_R1_e[nStripes_R1][nRadii]{};
0158   double x3a_R1[nStripes_R1][nRadii]{};
0159   double y3a_R1[nStripes_R1][nRadii]{};
0160   double x3a_R2[nStripes_R2][nRadii]{};
0161   double y3a_R2[nStripes_R2][nRadii]{};
0162   double x3a_R3[nStripes_R3][nRadii]{};
0163   double y3a_R3[nStripes_R3][nRadii]{};
0164
0165   // right midpoint - 3b
0166   double x3b_R1_e[nStripes_R1][nRadii]{};
0167   double y3b_R1_e[nStripes_R1][nRadii]{};
0168   double x3b_R1[nStripes_R1][nRadii]{};
0169   double y3b_R1[nStripes_R1][nRadii]{};
0170   double x3b_R2[nStripes_R2][nRadii]{};
0171   double y3b_R2[nStripes_R2][nRadii]{};
0172   double x3b_R3[nStripes_R3][nRadii]{};
0173   double y3b_R3[nStripes_R3][nRadii]{};
0174
0175   // Check which stripes get removed
0176   std::array<int, nRadii> nGoodStripes_R1_e = {};
0177   std::array<int, nRadii> nGoodStripes_R1 = {};
0178   std::array<int, nRadii> nGoodStripes_R2 = {};
0179   std::array<int, nRadii> nGoodStripes_R3 = {};
0180
0181   /// min stripe index
0182   static constexpr std::array<int, nRadii> keepThisAndAfter = {{1, 0, 1, 0, 1, 0, 1, 0}};
0183
0184   /// max stripe index
0185   static constexpr std::array<int, nRadii> keepUntil_R1_e = {{4, 4, 5, 4, 5, 5, 5, 5}};
0186   static constexpr std::array<int, nRadii> keepUntil_R1 = {{5, 5, 6, 5, 6, 5, 6, 5}};
0187   static constexpr std::array<int, nRadii> keepUntil_R2 = {{7, 7, 8, 7, 8, 8, 8, 8}};
0188   static constexpr std::array<int, nRadii> keepUntil_R3 = {{11, 10, 11, 11, 11, 11, 12, 11}};
0189
0190   std::array<int, nRadii> nStripesIn_R1_e = {};
0191   std::array<int, nRadii> nStripesIn_R1 = {};
0192   std::array<int, nRadii> nStripesIn_R2 = {};
0193   std::array<int, nRadii> nStripesIn_R3 = {};
0194   std::array<int, nRadii> nStripesBefore_R1_e = {};
0195   std::array<int, nRadii> nStripesBefore_R1 = {};
0196   std::array<int, nRadii> nStripesBefore_R2 = {};
0197   std::array<int, nRadii> nStripesBefore_R3 = {};
0198
0199   static constexpr int nStripesPerPetal = 213;
0200   static constexpr int nPetals = 18;
0201   static constexpr int nTotStripes = nStripesPerPetal * nPetals;
0202
0203   /// mean number of electrons per stripe
0204   int electrons_per_stripe = 300;
0205
0206   // number of electrons per deposited GeV in TPC gas
0207   /**
0208    * it is used to convert a given number of electrons into an energy
0209    * as expected by G4Hit. The energy is then converted back to a number of electrons
0210    * inside PHG4TpcElectronDrift
0211    */
0212   double electrons_per_gev = std::numeric_limits<double>::signaling_NaN();
0213
0214   /// delay between central membrane hits and trigger time (ns)
0215   int m_centralMembraneDelay = 0;
0216
0217   void CalculateVertices(
0218       int nStripes, int nPads,
0219       const std::array<double, nRadii>& R,
0220       std::array<double, nRadii>& spacing,
0221       double x1a[][nRadii], double y1a[][nRadii],
0222       double x1b[][nRadii], double y1b[][nRadii],
0223       double x2a[][nRadii], double y2a[][nRadii],
0224       double x2b[][nRadii], double y2b[][nRadii],
0225       double x3a[][nRadii], double y3a[][nRadii],
0226       double x3b[][nRadii], double y3b[][nRadii],
0227       double padfrac,
0228       double str_width[][nRadii],
0229       const std::array<double, nRadii>& widthmod,
0230       std::array<int, nRadii>& nGoodStripes,
0231       const std::array<int, nRadii>& keepUntil,
0232       std::array<int, nRadii>& nStripesIn,
0233       std::array<int, nRadii>& nStripesBefore);
0234
0235   static int SearchModule(int nStripes,
0236                           const double x1a[][nRadii], const double x1b[][nRadii],
0237                           const double x2a[][nRadii], const double x2b[][nRadii],
0238                           const double y1a[][nRadii], const double y1b[][nRadii],
0239                           const double y2a[][nRadii], const double y2b[][nRadii],
0240                           const double x3a[][nRadii], const double y3a[][nRadii],
0241                           const double x3b[][nRadii], const double y3b[][nRadii],
0242                           double x, double y, const std::array<int, nRadii>& nGoodStripes);
0243
0244   PHG4Hit* GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const;
0245 };
0246
0247 #endif