g4simulation/g4tpc/PHG4TpcCentralMembrane.cc

0001 #include "PHG4TpcCentralMembrane.h"
0002
0003 #include <phparameter/PHParameterInterface.h>  // for PHParameterInterface
0004
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4HitDefs.h>  // for get_volume_id
0008 #include <g4main/PHG4Hitv1.h>
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/SubsysReco.h>  // for SubsysReco
0012
0013 #include <phool/getClass.h>
0014 #include <phool/phool.h>  // for PHWHERE
0015
0016 #include <TVector3.h>
0017
0018 #include <iostream>  // for operator<<, endl, basi...
0019
0020 class PHCompositeNode;
0021
0022 namespace
0023 {
0024   // unique detector id for all direct lasers
0025   const int detId = PHG4HitDefs::get_volume_id("PHG4TpcCentralMembrane");
0026
0027   template <class T>
0028   constexpr T square(const T& x)
0029   {
0030     return x * x;
0031   }
0032
0033   template <class T>
0034   inline T get_r(const T& x, const T& y)
0035   {
0036     return std::sqrt(square(x) + square(y));
0037   }
0038
0039 }  // namespace
0040
0041 //_____________________________________________________________
0042 // all distances in mm, all angles in rad
0043 // class that generates stripes and dummy hit coordinates
0044 // stripes have width of one mm, length of one pad width, and are centered in middle of sector gaps
0045 PHG4TpcCentralMembrane::PHG4TpcCentralMembrane(const std::string& name)
0046   : SubsysReco(name)
0047   , PHParameterInterface(name)
0048 {
0049   InitializeParameters();
0050
0051   // set to 1.0 mm for all else
0052   for (int j = 0; j < nRadii; j++)
0053   {
0054     for (int i = 0; i < nStripes_R1; i++)
0055     {
0056       str_width_R1_e[i][j] = 1.0 * mm;
0057       str_width_R1[i][j] = 1.0 * mm;
0058     }
0059     for (auto& i : str_width_R2)
0060     {
0061       i[j] = 1.0 * mm;
0062     }
0063     for (auto& i : str_width_R3)
0064     {
0065       i[j] = 1.0 * mm;
0066     }
0067   }
0068
0069   CalculateVertices(nStripes_R1, nPads_R1, R1_e, spacing_R1_e, x1a_R1_e, y1a_R1_e, x1b_R1_e, y1b_R1_e, x2a_R1_e, y2a_R1_e, x2b_R1_e, y2b_R1_e, x3a_R1_e, y3a_R1_e, x3b_R1_e, y3b_R1_e, padfrac_R1, str_width_R1_e, widthmod_R1_e, nGoodStripes_R1_e, keepUntil_R1_e, nStripesIn_R1_e, nStripesBefore_R1_e);
0070   CalculateVertices(nStripes_R1, nPads_R1, R1, spacing_R1, x1a_R1, y1a_R1, x1b_R1, y1b_R1, x2a_R1, y2a_R1, x2b_R1, y2b_R1, x3a_R1, y3a_R1, x3b_R1, y3b_R1, padfrac_R1, str_width_R1, widthmod_R1, nGoodStripes_R1, keepUntil_R1, nStripesIn_R1, nStripesBefore_R1);
0071   CalculateVertices(nStripes_R2, nPads_R2, R2, spacing_R2, x1a_R2, y1a_R2, x1b_R2, y1b_R2, x2a_R2, y2a_R2, x2b_R2, y2b_R2, x3a_R2, y3a_R2, x3b_R2, y3b_R2, padfrac_R2, str_width_R2, widthmod_R2, nGoodStripes_R2, keepUntil_R2, nStripesIn_R2, nStripesBefore_R2);
0072   CalculateVertices(nStripes_R3, nPads_R3, R3, spacing_R3, x1a_R3, y1a_R3, x1b_R3, y1b_R3, x2a_R3, y2a_R3, x2b_R3, y2b_R3, x3a_R3, y3a_R3, x3b_R3, y3b_R3, padfrac_R3, str_width_R3, widthmod_R3, nGoodStripes_R3, keepUntil_R3, nStripesIn_R3, nStripesBefore_R3);
0073 }
0074
0075 //______________________________________________________
0076 PHG4TpcCentralMembrane::~PHG4TpcCentralMembrane()
0077 {
0078   for (auto&& hit : PHG4Hits)
0079   {
0080     delete hit;
0081   }
0082 }
0083
0084 //_____________________________________________________________
0085 int PHG4TpcCentralMembrane::InitRun(PHCompositeNode* /* topNode */)
0086 {
0087   // setup parameters
0088   UpdateParametersWithMacro();
0089   electrons_per_stripe = get_int_param("electrons_per_stripe");
0090   electrons_per_gev = get_double_param("electrons_per_gev");
0091
0092   std::cout << "PHG4TpcCentralMembrane::InitRun - electrons_per_stripe: " << electrons_per_stripe << std::endl;
0093   std::cout << "PHG4TpcCentralMembrane::InitRun - electrons_per_gev " << electrons_per_gev << std::endl;
0094
0095   // reset g4hits
0096   for (auto&& hit : PHG4Hits)
0097   {
0098     delete hit;
0099   }
0100
0101   PHG4Hits.clear();
0102
0103   /*
0104    * utility function to
0105    * - duplicate generated G4Hit to cover both sides of the central membrane
0106    * - adjust hit time and z,
0107    * - insert in container
0108    */
0109   auto adjust_hits = [&](PHG4Hit* source)
0110   {
0111     // adjust time to account for central membrane delay
0112     source->set_t(0, m_centralMembraneDelay);
0113     source->set_t(1, m_centralMembraneDelay);
0114
0115     // assign to positive side
0116     source->set_z(0, 1.);
0117     source->set_z(1, 1.);
0118     PHG4Hits.push_back(source);
0119
0120     // clone
0121     // assign to negative side and insert in list
0122     auto* copy = new PHG4Hitv1(source);
0123     copy->set_z(0, -1.);
0124     copy->set_z(1, -1.);
0125     PHG4Hits.push_back(copy);
0126   };
0127
0128   // loop over petalID
0129   for (int i = 0; i < 18; i++)
0130   {
0131     // loop over radiusID
0132     for (int j = 0; j < 8; j++)
0133     {
0134       // loop over stripeID
0135       for (int k = 0; k < nGoodStripes_R1_e[j]; k++)
0136       {
0137         adjust_hits(GetPHG4HitFromStripe(i, 0, j, k, electrons_per_stripe));
0138       }
0139
0140       // loop over stripeID
0141       for (int k = 0; k < nGoodStripes_R1[j]; k++)
0142       {
0143         adjust_hits(GetPHG4HitFromStripe(i, 1, j, k, electrons_per_stripe));
0144       }
0145
0146       // loop over stripeID
0147       for (int k = 0; k < nGoodStripes_R2[j]; k++)
0148       {
0149         adjust_hits(GetPHG4HitFromStripe(i, 2, j, k, electrons_per_stripe));
0150       }
0151
0152       // loop over stripeID
0153       for (int k = 0; k < nGoodStripes_R3[j]; k++)
0154       {
0155         adjust_hits(GetPHG4HitFromStripe(i, 3, j, k, electrons_per_stripe));
0156       }
0157     }
0158   }
0159
0160   m_eventNum = 0;
0161
0162   return Fun4AllReturnCodes::EVENT_OK;
0163 }
0164
0165 //_____________________________________________________________
0166 int PHG4TpcCentralMembrane::process_event(PHCompositeNode* topNode)
0167 {
0168   if (m_eventNum % m_eventModulo != 0)
0169   {
0170     if (Verbosity())
0171     {
0172       std::cout << "Event " << m_eventNum << " will not generate CM hits" << std::endl;
0173     }
0174     m_eventNum++;
0175     return Fun4AllReturnCodes::EVENT_OK;
0176   }
0177
0178   // load g4hit container
0179   auto* g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0180   if (!g4hitcontainer)
0181   {
0182     std::cout << PHWHERE << "Could not locate g4 hit node " << hitnodename << std::endl;
0183     return Fun4AllReturnCodes::ABORTRUN;
0184   }
0185
0186   // copy all hits from G4hits vector into container
0187   for (const auto& hit : PHG4Hits)
0188   {
0189     auto* copy = new PHG4Hitv1(hit);
0190     g4hitcontainer->AddHit(detId, copy);
0191   }
0192
0193   m_eventNum++;
0194
0195   return Fun4AllReturnCodes::EVENT_OK;
0196 }
0197
0198 //_____________________________________________________________
0199 void PHG4TpcCentralMembrane::SetDefaultParameters()
0200 {
0201   // same gas parameters as in PHG4TpcElectronDrift::SetDefaultParameters
0202
0203   // Data on gasses @20 C and 760 Torr from the following source:
0204   // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
0205   // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
0206   // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
0207   static constexpr double Ne_dEdx = 1.56;    // keV/cm
0208   static constexpr double CF4_dEdx = 7.00;   // keV/cm
0209   static constexpr double Ne_NTotal = 43;    // Number/cm
0210   static constexpr double CF4_NTotal = 100;  // Number/cm
0211   static constexpr double Tpc_NTot = (0.5 * Ne_NTotal) + (0.5 * CF4_NTotal);
0212   static constexpr double Tpc_dEdx = (0.5 * Ne_dEdx) + (0.5 * CF4_dEdx);
0213   static constexpr double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
0214
0215   // number of electrons per deposited GeV in TPC gas
0216   set_default_double_param("electrons_per_gev", Tpc_ElectronsPerKeV * 1000000.);
0217
0218   /// mean number of electrons per stripe
0219   set_default_int_param("electrons_per_stripe", 100);
0220 }
0221
0222 //_____________________________________________________________
0223 void PHG4TpcCentralMembrane::CalculateVertices(
0224     int nStripes, int nPads,
0225     const std::array<double, nRadii>& R,
0226     std::array<double, nRadii>& spacing,
0227     double x1a[][nRadii], double y1a[][nRadii],
0228     double x1b[][nRadii], double y1b[][nRadii],
0229     double x2a[][nRadii], double y2a[][nRadii],
0230     double x2b[][nRadii], double y2b[][nRadii],
0231     double x3a[][nRadii], double y3a[][nRadii],
0232     double x3b[][nRadii], double y3b[][nRadii],
0233     double padfrac,
0234     double str_width[][nRadii],
0235     const std::array<double, nRadii>& widthmod,
0236     std::array<int, nRadii>& nGoodStripes,
0237     const std::array<int, nRadii>& keepUntil,
0238     std::array<int, nRadii>& nStripesIn,
0239     std::array<int, nRadii>& nStripesBefore)
0240 {
0241   const double phi_module = M_PI / 6.0;  // angle span of a module
0242   const int pr_mult = 3;                 // multiples of intrinsic resolution of pads
0243   const int dw_mult = 8;                 // multiples of diffusion width
0244   const double diffwidth = 0.6 * mm;     // diffusion width
0245   const double adjust = 0.015;           // arbitrary angle to center the pattern in a petal
0246
0247   double theta = 0.0;
0248   // center coords
0249   std::vector<std::vector<double>> cx(nStripes, std::vector<double>(nRadii));
0250   std::vector<std::vector<double>> cy(nStripes, std::vector<double>(nRadii));
0251   // corner coords
0252   /* double tempX1a[nStripes][nRadii], tempY1a[nStripes][nRadii];
0253   double tempX1b[nStripes][nRadii], tempY1b[nStripes][nRadii];
0254   double tempX2a[nStripes][nRadii], tempY2a[nStripes][nRadii];
0255   double tempX2b[nStripes][nRadii], tempY2b[nStripes][nRadii];
0256   double rotatedX1a[nStripes][nRadii], rotatedY1a[nStripes][nRadii];
0257   double rotatedX1b[nStripes][nRadii], rotatedY1b[nStripes][nRadii];
0258   double rotatedX2a[nStripes][nRadii], rotatedY2a[nStripes][nRadii];
0259   double rotatedX2b[nStripes][nRadii], rotatedY2b[nStripes][nRadii]; */
0260
0261   // calculate spacing first:
0262   for (int i = 0; i < nRadii; i++)
0263   {
0264     spacing[i] = 2.0 * ((dw_mult * diffwidth / R[i]) + (pr_mult * phi_module / nPads));
0265   }
0266
0267   // vertex calculation
0268   for (int j = 0; j < nRadii; j++)
0269   {
0270     int i_out = 0;
0271     for (int i = keepThisAndAfter[j]; i < keepUntil[j]; i++)
0272     {
0273       if (j % 2 == 0)
0274       {
0275         theta = i * spacing[j] + (spacing[j] / 2) - adjust;
0276         cx[i_out][j] = R[j] * cos(theta);
0277         cy[i_out][j] = R[j] * sin(theta);
0278       }
0279       else
0280       {
0281         theta = (i + 1) * spacing[j] - adjust;
0282         cx[i_out][j] = R[j] * cos(theta);
0283         cy[i_out][j] = R[j] * sin(theta);
0284       }
0285
0286       TVector3 corner[4];
0287       corner[0].SetXYZ(-padfrac + arc_r, -(widthmod[j] * str_width[i][j]) / 2, 0);  //"1a" = length of the pad, but not including the arc piece
0288       corner[1].SetXYZ(padfrac - arc_r, -(widthmod[j] * str_width[i][j]) / 2, 0);   //"1b" = length of the pad, but not including the arc piece
0289       corner[2].SetXYZ(-padfrac + arc_r, (widthmod[j] * str_width[i][j]) / 2, 0);   //"2a" = length of the pad, but not including the arc piece
0290       corner[3].SetXYZ(padfrac - arc_r, (widthmod[j] * str_width[i][j]) / 2, 0);    //"2b" = length of the pad, but not including the arc piece
0291
0292       TVector3 rotatedcorner[4];
0293       for (int n = 0; n < 4; n++)
0294       {
0295         rotatedcorner[n] = corner[n];
0296         rotatedcorner[n].RotateZ(theta);
0297       }
0298
0299       x1a[i_out][j] = rotatedcorner[0].X() + cx[i_out][j];
0300       x1b[i_out][j] = rotatedcorner[1].X() + cx[i_out][j];
0301       x2a[i_out][j] = rotatedcorner[2].X() + cx[i_out][j];
0302       x2b[i_out][j] = rotatedcorner[3].X() + cx[i_out][j];
0303
0304       y1a[i_out][j] = rotatedcorner[0].Y() + cy[i_out][j];
0305       y1b[i_out][j] = rotatedcorner[1].Y() + cy[i_out][j];
0306       y2a[i_out][j] = rotatedcorner[2].Y() + cy[i_out][j];
0307       y2b[i_out][j] = rotatedcorner[3].Y() + cy[i_out][j];
0308
0309       /* x1a[i_out][j] = cx[i_out][j] - padfrac + arc_r;
0310       y1a[i_out][j] = cy[i_out][j] - str_width/2;
0311       x1b[i_out][j] = cx[i_out][j] + padfrac - arc_r;
0312       y1b[i_out][j] = cy[i_out][j] - str_width/2;
0313       x2a[i_out][j] = cx[i_out][j] - padfrac + arc_r;
0314       y2a[i_out][j] = cy[i_out][j] + str_width/2;
0315       x2b[i_out][j] = cx[i_out][j] + padfrac - arc_r;
0316       y2b[i_out][j] = cy[i_out][j] + str_width/2;
0317
0318       tempX1a[i_out][j] = x1a[i_out][j] - cx[i_out][j];
0319       tempY1a[i_out][j] = y1a[i_out][j] - cy[i_out][j];
0320       tempX1b[i_out][j] = x1b[i_out][j] - cx[i_out][j];
0321       tempY1b[i_out][j] = y1b[i_out][j] - cy[i_out][j];
0322       tempX2a[i_out][j] = x2a[i_out][j] - cx[i_out][j];
0323       tempY2a[i_out][j] = y2a[i_out][j] - cy[i_out][j];
0324       tempX2b[i_out][j] = x2b[i_out][j] - cx[i_out][j];
0325       tempY2b[i_out][j] = y2b[i_out][j] - cy[i_out][j];
0326
0327       rotatedX1a[i_out][j] = tempX1a[i_out][j]*cos(theta) - tempY1a[i_out][j]*sin(theta);
0328       rotatedY1a[i_out][j] = tempX1a[i_out][j]*sin(theta) + tempY1a[i_out][j]*cos(theta);
0329       rotatedX1b[i_out][j] = tempX1b[i_out][j]*cos(theta) - tempY1b[i_out][j]*sin(theta);
0330       rotatedY1b[i_out][j] = tempX1b[i_out][j]*sin(theta) + tempY1b[i_out][j]*cos(theta);
0331       rotatedX2a[i_out][j] = tempX2a[i_out][j]*cos(theta) - tempY2a[i_out][j]*sin(theta);
0332       rotatedY2a[i_out][j] = tempX2a[i_out][j]*sin(theta) + tempY2a[i_out][j]*cos(theta);
0333       rotatedX2b[i_out][j] = tempX2b[i_out][j]*cos(theta) - tempY2b[i_out][j]*sin(theta);
0334       rotatedY2b[i_out][j] = tempX2b[i_out][j]*sin(theta) + tempY2b[i_out][j]*cos(theta);*/
0335
0336       /* x1a[i_out][j] = rotatedX1a[i_out][j] + cx[i_out][j];
0337       y1a[i_out][j] = rotatedY1a[i_out][j] + cy[i_out][j];
0338       x1b[i_out][j] = rotatedX1b[i_out][j] + cx[i_out][j];
0339       y1b[i_out][j] = rotatedY1b[i_out][j] + cy[i_out][j];
0340       x2a[i_out][j] = rotatedX2a[i_out][j] + cx[i_out][j];
0341       y2a[i_out][j] = rotatedY2a[i_out][j] + cy[i_out][j];
0342       x2b[i_out][j] = rotatedX2b[i_out][j] + cx[i_out][j];
0343       y2b[i_out][j] = rotatedY2b[i_out][j] + cy[i_out][j]; */
0344
0345       x3a[i_out][j] = (x1a[i_out][j] + x2a[i_out][j]) / 2.0;
0346       y3a[i_out][j] = (y1a[i_out][j] + y2a[i_out][j]) / 2.0;
0347       x3b[i_out][j] = (x1b[i_out][j] + x2b[i_out][j]) / 2.0;
0348       y3b[i_out][j] = (y1b[i_out][j] + y2b[i_out][j]) / 2.0;
0349
0350       i_out++;
0351
0352       nStripesBefore_R1_e[0] = 0;
0353
0354       nStripesIn[j] = keepUntil[j] - keepThisAndAfter[j];
0355       if (j == 0)
0356       {
0357         nStripesBefore[j] = 0;
0358       }
0359       else
0360       {
0361         nStripesBefore[j] = nStripesIn[j - 1] + nStripesBefore[j - 1];
0362       }
0363       nStripesBefore_R1_e[0] = 0;
0364     }
0365     nGoodStripes[j] = i_out;
0366   }
0367 }
0368
0369 int PHG4TpcCentralMembrane::SearchModule(int /*nStripes*/,
0370                                          const double x1a[][nRadii], const double x1b[][nRadii],
0371                                          const double x2a[][nRadii], const double x2b[][nRadii],
0372                                          const double y1a[][nRadii], const double y1b[][nRadii],
0373                                          const double y2a[][nRadii], const double y2b[][nRadii],
0374                                          const double x3a[][nRadii], const double y3a[][nRadii],
0375                                          const double x3b[][nRadii], const double y3b[][nRadii],
0376                                          double x, double y, const std::array<int, nRadii>& nGoodStripes)
0377 {
0378   int c = 0;
0379
0380   for (int j = 0; j < nRadii; j++)
0381   {
0382     for (int i = 0; i < nGoodStripes[j]; i++)
0383     {
0384       if (((y1a[i][j] > y) != (y2a[i][j] > y) && (x < (x2a[i][j] - x1a[i][j]) * (y - y1a[i][j]) / (y2a[i][j] - y1a[i][j]) + x1a[i][j])))
0385       {
0386         c = !c;
0387       }
0388       if (((y1b[i][j] > y) != (y1a[i][j] > y) && (x < (x1a[i][j] - x1b[i][j]) * (y - y1b[i][j]) / (y1a[i][j] - y1b[i][j]) + x1b[i][j])))
0389       {
0390         c = !c;
0391       }
0392       if (((y2b[i][j] > y) != (y1b[i][j] > y) && (x < (x1b[i][j] - x2b[i][j]) * (y - y2b[i][j]) / (y1b[i][j] - y2b[i][j]) + x2b[i][j])))
0393       {
0394         c = !c;
0395       }
0396       if (((y2a[i][j] > y) != (y2b[i][j] > y) && (x < (x2b[i][j] - x2a[i][j]) * (y - y2a[i][j]) / (y2b[i][j] - y2a[i][j]) + x2a[i][j])))
0397       {
0398         c = !c;
0399       }
0400
0401       // check inside arcs
0402       if (c == 0)
0403       {
0404         if (((x - x3a[i][j]) * (x - x3a[i][j]) + (y - y3a[i][j]) * (y - y3a[i][j])) <= arc_r * arc_r || ((x - x3b[i][j]) * (x - x3b[i][j]) + (y - y3b[i][j]) * (y - y3b[i][j])) <= arc_r * arc_r)
0405         {
0406           c = !c;
0407         }
0408       }
0409     }
0410   }
0411   return c;
0412 }
0413
0414 int PHG4TpcCentralMembrane::getSearchResult(double xcheck, double ycheck) const
0415 {
0416   const double phi_petal = M_PI / 9.0;  // angle span of one petal
0417   const double end_R1_e = 312.0 * mm;   // arbitrary radius between R1_e and R1
0418   const double end_R1 = 408.0 * mm;     // arbitrary radius between R1 and R2
0419   const double end_R2 = 580.0 * mm;     // arbitrary radius between R2 and R3
0420
0421   double r;
0422   double phi;
0423   double phimod;
0424   double xmod;
0425   double ymod;
0426
0427   r = sqrt((xcheck * xcheck) + (ycheck * ycheck));
0428   phi = atan(ycheck / xcheck);
0429   if ((xcheck < 0.0) && (ycheck > 0.0))
0430   {
0431     phi = phi + M_PI;
0432   }
0433   else if ((xcheck > 0.0) && (ycheck < 0.0))
0434   {
0435     phi = phi + 2.0 * M_PI;
0436   }
0437
0438   phimod = fmod(phi, phi_petal);
0439   xmod = r * cos(phimod);
0440   ymod = r * sin(phimod);
0441
0442   int result = 0;
0443
0444   if (r <= end_R1_e)
0445   {
0446     result = SearchModule(nStripes_R1, x1a_R1_e, x1b_R1_e, x2a_R1_e, x2b_R1_e, y1a_R1_e, y1b_R1_e, y2a_R1_e, y2b_R1_e, x3a_R1_e, y3a_R1_e, x3b_R1_e, y3b_R1_e, xmod, ymod, nGoodStripes_R1_e);
0447   }
0448   else if ((r > end_R1_e) && (r <= end_R1))
0449   {
0450     result = SearchModule(nStripes_R1, x1a_R1, x1b_R1, x2a_R1, x2b_R1, y1a_R1, y1b_R1, y2a_R1, y2b_R1, x3a_R1, y3a_R1, x3b_R1, y3b_R1, xmod, ymod, nGoodStripes_R1);
0451   }
0452   else if ((r > end_R1) && (r <= end_R2))
0453   {
0454     result = SearchModule(nStripes_R2, x1a_R2, x1b_R2, x2a_R2, x2b_R2, y1a_R2, y1b_R2, y2a_R2, y2b_R2, x3a_R2, y3a_R2, x3b_R2, y3b_R2, xmod, ymod, nGoodStripes_R2);
0455   }
0456   else if ((r > end_R2) && (r <= end_CM))
0457   {
0458     result = SearchModule(nStripes_R3, x1a_R3, x1b_R3, x2a_R3, x2b_R3, y1a_R3, y1b_R3, y2a_R3, y2b_R3, x3a_R3, y3a_R3, x3b_R3, y3b_R3, xmod, ymod, nGoodStripes_R3);
0459   }
0460
0461   return result;
0462 }
0463
0464 PHG4Hit* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const
0465 {                                       // this function generates a PHG4 hit using coordinates from a stripe
0466   const double phi_petal = M_PI / 9.0;  // angle span of one petal
0467   PHG4Hit* hit;
0468   TVector3 dummyPos0;
0469   TVector3 dummyPos1;
0470
0471   // could put in some sanity checks here but probably not necessary since this is only really used within the class
0472   // petalID ranges 0-17, module ID 0-3, stripeID varies - nGoodStripes for each module
0473   // radiusID ranges 0-7
0474
0475   // from phg4tpcsteppingaction.cc
0476   hit = new PHG4Hitv1();
0477   hit->set_layer(-1);  // dummy number
0478   // here we set the entrance values in cm
0479   if (moduleID == 0)
0480   {
0481     hit->set_x(0, x3a_R1_e[stripeID][radiusID] / cm);
0482     hit->set_y(0, y3a_R1_e[stripeID][radiusID] / cm);
0483   }
0484   else if (moduleID == 1)
0485   {
0486     hit->set_x(0, x3a_R1[stripeID][radiusID] / cm);
0487     hit->set_y(0, y3a_R1[stripeID][radiusID] / cm);
0488   }
0489   else if (moduleID == 2)
0490   {
0491     hit->set_x(0, x3a_R2[stripeID][radiusID] / cm);
0492     hit->set_y(0, y3a_R2[stripeID][radiusID] / cm);
0493   }
0494   else if (moduleID == 3)
0495   {
0496     hit->set_x(0, x3a_R3[stripeID][radiusID] / cm);
0497     hit->set_y(0, y3a_R3[stripeID][radiusID] / cm);
0498   }
0499   hit->set_z(0, 0.0 / cm);
0500
0501   // check if you need to rotate coords to another petal
0502   if (petalID > 0)
0503   {
0504     dummyPos0.SetXYZ(hit->get_x(0), hit->get_y(0), hit->get_z(0));
0505     dummyPos0.RotateZ(petalID * phi_petal);
0506     hit->set_x(0, dummyPos0.X());
0507     hit->set_y(0, dummyPos0.Y());
0508   }
0509
0510   // TODO: use actual stripe direction for the momentum
0511   hit->set_px(1, 500.0);
0512   hit->set_py(1, 500.0);
0513   hit->set_pz(1, 500.0);
0514
0515   // time in ns
0516   hit->set_t(0, 0);
0517
0518   // set and save the track ID
0519   hit->set_trkid(-1);  // dummy number
0520
0521   // here we just update the exit values, it will be overwritten
0522   // for every step until we leave the volume or the particle
0523   // ceases to exist
0524   if (moduleID == 0)
0525   {
0526     hit->set_x(1, x3b_R1_e[stripeID][radiusID] / cm);
0527     hit->set_y(1, y3b_R1_e[stripeID][radiusID] / cm);
0528   }
0529   else if (moduleID == 1)
0530   {
0531     hit->set_x(1, x3b_R1[stripeID][radiusID] / cm);
0532     hit->set_y(1, y3b_R1[stripeID][radiusID] / cm);
0533   }
0534   else if (moduleID == 2)
0535   {
0536     hit->set_x(1, x3b_R2[stripeID][radiusID] / cm);
0537     hit->set_y(1, y3b_R2[stripeID][radiusID] / cm);
0538   }
0539   else if (moduleID == 3)
0540   {
0541     hit->set_x(1, x3b_R3[stripeID][radiusID] / cm);
0542     hit->set_y(1, y3b_R3[stripeID][radiusID] / cm);
0543   }
0544   hit->set_z(1, 0.0 / cm);
0545
0546   // check if you need to rotate coords to another petal
0547   if (petalID > 0)
0548   {
0549     dummyPos1.SetXYZ(hit->get_x(1), hit->get_y(1), hit->get_z(1));
0550     dummyPos1.RotateZ(petalID * phi_petal);
0551     hit->set_x(1, dummyPos1.X());
0552     hit->set_y(1, dummyPos1.Y());
0553   }
0554
0555   // TODO: use actual stripe direction for the momentum
0556   hit->set_px(1, 500.0);
0557   hit->set_py(1, 500.0);
0558   hit->set_pz(1, 500.0);
0559
0560   hit->set_t(1, 0);  // dummy number, nanosecond
0561
0562   // calculate deposited energy corresponding to number of electrons per stripe
0563   const double edep = nElectrons / electrons_per_gev;
0564   hit->set_edep(edep);
0565   hit->set_eion(edep);
0566
0567   /*
0568   if (hit->get_edep()){ //print out hits
0569     double rin = sqrt(hit->get_x(0) * hit->get_x(0) + hit->get_y(0) * hit->get_y(0));
0570     double rout = sqrt(hit->get_x(1) * hit->get_x(1) + hit->get_y(1) * hit->get_y(1));
0571     std::cout << "Added Tpc g4hit with rin, rout = " << rin << "  " << rout
0572          << " g4hitid " << hit->get_hit_id() << std::endl;
0573     std::cout << " xin " << hit->get_x(0)
0574          << " yin " << hit->get_y(0)
0575          << " zin " << hit->get_z(0)
0576          << " rin " << rin
0577          << std::endl;
0578     std::cout << " xout " << hit->get_x(1)
0579          << " yout " << hit->get_y(1)
0580          << " zout " << hit->get_z(1)
0581          << " rout " << rout
0582          << std::endl;
0583     std::cout << " xav " << (hit->get_x(1) + hit->get_x(0)) / 2.0
0584          << " yav " << (hit->get_y(1) + hit->get_y(0)) / 2.0
0585          << " zav " << (hit->get_z(1) + hit->get_z(0)) / 2.0
0586          << " rav " << (rout + rin) / 2.0
0587          << std::endl;
0588   }
0589   */
0590
0591   return hit;
0592 }
0593
0594 int PHG4TpcCentralMembrane::getStripeID(double xcheck, double ycheck) const
0595 {
0596   // check if point came from stripe then see which stripe it is
0597   // 213 stripes in a petal, 18 petals, ntotstripes = 3834
0598   int result;
0599   int fullID = -1;
0600   // const double adjust = 0.015; //arbitrary angle to center the pattern in a petal
0601   const double phi_petal = M_PI / 9.0;  // angle span of one petal
0602
0603   // check if in a stripe
0604   result = getSearchResult(xcheck, ycheck);
0605
0606   // find which stripe
0607   if (result == 1)
0608   {
0609     // std::cout << "on a stripe" << std::endl;
0610     // convert coords to radius n angle
0611     double r = sqrt((xcheck * xcheck) + (ycheck * ycheck));
0612     double phi = atan(ycheck / xcheck);
0613     if ((xcheck < 0.0) && (ycheck > 0.0))
0614     {
0615       phi = phi + M_PI;
0616     }
0617     else if ((xcheck > 0.0) && (ycheck < 0.0))
0618     {
0619       phi = phi + 2.0 * M_PI;
0620     }
0621     // get angle within first petal
0622     double phimod = fmod(phi, phi_petal);
0623     double xmod = r * cos(phimod);
0624     double ymod = r * sin(phimod);
0625
0626     int petalID = phi / phi_petal;
0627
0628     int phiID = 0;
0629     for (int j = 0; j < nRadii; j++)
0630     {
0631       if (((R1_e[j] - padfrac_R1) < r) && (r < (R1_e[j] + padfrac_R1)))
0632       {  // check if radius is in stripe
0633         int rID = j;
0634         std::cout << "rID: " << rID << std::endl;
0635         //'angle' is to the center of a stripe
0636         for (int i = 0; i < nGoodStripes_R1_e[j]; i++)
0637         {
0638           // if (j % 2 == 0){
0639           // theta = i*spacing[j];
0640           // angle = theta + (spacing[j]/2) - adjust;
0641           //  look at distance from center line of stripe
0642           //  if distance from x,y to center line < str_width
0643           //  calculate slope n then do dist
0644
0645           double m = (y3b_R1_e[i][j] - y3a_R1_e[i][j]) / (x3b_R1_e[i][j] - x3a_R1_e[i][j]);
0646           /*std::cout << "y2: " << y3b_R1_e[i][j] << std::endl;
0647     std::cout << "y1: " << y3a_R1_e[i][j] << std::endl;
0648     std::cout << "x2: " << x3b_R1_e[i][j] << std::endl;
0649     std::cout << "x1: " << x3a_R1_e[i][j] << std::endl;
0650     std::cout << "xc: " << xcheck << std::endl;
0651     std::cout << "yc: " << ycheck << std::endl;
0652           std::cout << "m: " << m << std::endl;  */
0653           // std::cout << fabs((-m)*xcheck + ycheck) << std::endl;
0654           double dist = fabs(((-m) * xmod) + ymod) / sqrt(1 + (m * m));
0655           // std::cout << "dist:" << dist << std::endl;
0656           if (dist < ((widthmod_R1_e[j] * str_width_R1_e[i][j]) / 2.0))
0657           {
0658             phiID = i;
0659             // std::cout << "phiID: " << phiID << std::endl;
0660           }
0661         }
0662
0663         std::cout << "nStripesBefore: " << nStripesBefore_R1_e[j] << std::endl;
0664         fullID = petalID * nStripesPerPetal + nStripesBefore_R1_e[j] + phiID;
0665         // std::cout << "fullID: " << fullID << std::endl;
0666       }
0667       else if (((R1[j] - padfrac_R1) < r) && (r < (R1[j] + padfrac_R1)))
0668       {
0669         // std::cout << "R1" << std::endl;
0670         for (int i = 0; i < nGoodStripes_R1[j]; i++)
0671         {
0672           // look at distance from center line of stripe
0673           double m = (y3b_R1[i][j] - y3a_R1[i][j]) / (x3b_R1[i][j] - x3a_R1[i][j]);
0674           double dist = fabs((m * xmod) - ymod) / sqrt(1 + (m * m));
0675           if (dist < ((widthmod_R1[j] * str_width_R1[i][j]) / 2.0))
0676           {
0677             phiID = i;
0678           }
0679         }
0680
0681         fullID = petalID * nStripesPerPetal + nStripesBefore_R1[j] + phiID;
0682       }
0683       else if (((R2[j] - padfrac_R2) < r) && (r < (R2[j] + padfrac_R2)))
0684       {
0685         // std::cout << "R2" << std::endl;
0686         for (int i = 0; i < nGoodStripes_R2[j]; i++)
0687         {
0688           // look at distance from center line of stripe
0689           double m = (y3b_R2[i][j] - y3a_R2[i][j]) / (x3b_R2[i][j] - x3a_R2[i][j]);
0690           double dist = fabs((m * xmod) - ymod) / sqrt(1 + (m * m));
0691           if (dist < ((widthmod_R2[j] * str_width_R2[i][j]) / 2.0))
0692           {
0693             phiID = i;
0694           }
0695         }
0696
0697         fullID = petalID * nStripesPerPetal + nStripesBefore_R2[j] + phiID;
0698       }
0699       else if (((R3[j] - padfrac_R3) < r) && (r < (R3[j] + padfrac_R3)))
0700       {
0701         // std::cout << "R3" << std::endl;
0702         for (int i = 0; i < nGoodStripes_R3[j]; i++)
0703         {
0704           // look at distance from center line of stripe
0705           double m = (y3b_R3[i][j] - y3a_R3[i][j]) / (x3b_R3[i][j] - x3a_R3[i][j]);
0706           double dist = fabs((m * xmod) - ymod) / sqrt(1 + (m * m));
0707           if (dist < ((widthmod_R3[j] * str_width_R3[i][j]) / 2.0))
0708           {
0709             phiID = i;
0710           }
0711         }
0712
0713         fullID = petalID * nStripesPerPetal + nStripesBefore_R3[j] + phiID;
0714       }
0715     }
0716   }
0717   else
0718   {
0719     fullID = -1;
0720   }
0721
0722   return fullID;
0723 }