ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcCentralMembrane.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcCentralMembrane.cc
2 
3 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
4 
5 #include <g4main/PHG4Hit.h>
7 #include <g4main/PHG4HitDefs.h> // for get_volume_id
8 #include <g4main/PHG4Hitv1.h>
9 
11 #include <fun4all/SubsysReco.h> // for SubsysReco
12 
13 #include <phool/getClass.h>
14 #include <phool/phool.h> // for PHWHERE
15 
16 #include <TVector3.h>
17 
18 #include <iostream> // for operator<<, endl, basi...
19 
20 class PHCompositeNode;
21 
22 namespace
23 {
24  // unique detector id for all direct lasers
25  static const int detId = PHG4HitDefs::get_volume_id("PHG4TpcCentralMembrane");
26 
27 } // namespace
28 
29 //_____________________________________________________________
30 // all distances in mm, all angles in rad
31 // class that generates stripes and dummy hit coordinates
32 // stripes have width of one mm, length of one pad width, and are centered in middle of sector gaps
34  : SubsysReco(name)
35  , PHParameterInterface(name)
36 {
38 
39  // set to 1.0 mm for all else
40  for (int j = 0; j < nRadii; j++)
41  {
42  for (int i = 0; i < nStripes_R1; i++)
43  {
44  str_width_R1_e[i][j] = 1.0 * mm;
45  str_width_R1[i][j] = 1.0 * mm;
46  }
47  for (int i = 0; i < nStripes_R2; i++)
48  {
49  str_width_R2[i][j] = 1.0 * mm;
50  }
51  for (int i = 0; i < nStripes_R3; i++)
52  {
53  str_width_R3[i][j] = 1.0 * mm;
54  }
55  }
56 
61 }
62 
63 //______________________________________________________
65 {
66  for (auto&& hit : PHG4Hits)
67  {
68  delete hit;
69  }
70  for (auto&& hit : BotVertices)
71  {
72  delete hit;
73  }
74  for (auto&& hit : TopVertices)
75  {
76  delete hit;
77  }
78 }
79 
80 //_____________________________________________________________
82 {
83  // setup parameters
85  electrons_per_stripe = get_int_param("electrons_per_stripe");
86  electrons_per_gev = get_double_param("electrons_per_gev");
87 
88  std::cout << "PHG4TpcCentralMembrane::InitRun - electrons_per_stripe: " << electrons_per_stripe << std::endl;
89  std::cout << "PHG4TpcCentralMembrane::InitRun - electrons_per_gev " << electrons_per_gev << std::endl;
90 
91  // make sure G4Hit container exists
92  hitnodename = "G4HIT_" + detector;
93  auto* g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
94  if (!g4hit)
95  {
96  std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
98  }
99 
100  // reset vertices and g4hits
101  for (auto&& hit : PHG4Hits)
102  {
103  delete hit;
104  }
105  PHG4Hits.clear();
106 
107  for (auto&& hit : BotVertices)
108  {
109  delete hit;
110  }
111  BotVertices.clear();
112 
113  for (auto&& hit : TopVertices)
114  {
115  delete hit;
116  }
117  TopVertices.clear();
118 
119  for (int i = 0; i < 18; i++)
120  { // loop over petalID
121  for (int j = 0; j < 8; j++)
122  { // loop over radiusID
123  for (int k = 0; k < nGoodStripes_R1_e[j]; k++)
124  { // loop over stripeID
125  PHG4Hits.push_back(GetPHG4HitFromStripe(i, 0, j, k, electrons_per_stripe));
126  BotVertices.push_back(GetBotVerticesFromStripe(0, j, k));
127  TopVertices.push_back(GetTopVerticesFromStripe(0, j, k));
128  }
129  for (int k = 0; k < nGoodStripes_R1[j]; k++)
130  { // loop over stripeID
131  PHG4Hits.push_back(GetPHG4HitFromStripe(i, 1, j, k, electrons_per_stripe));
132  BotVertices.push_back(GetBotVerticesFromStripe(1, j, k));
133  TopVertices.push_back(GetTopVerticesFromStripe(1, j, k));
134  }
135  for (int k = 0; k < nGoodStripes_R2[j]; k++)
136  { // loop over stripeID
137  PHG4Hits.push_back(GetPHG4HitFromStripe(i, 2, j, k, electrons_per_stripe));
138  BotVertices.push_back(GetBotVerticesFromStripe(2, j, k));
139  TopVertices.push_back(GetTopVerticesFromStripe(2, j, k));
140  }
141  for (int k = 0; k < nGoodStripes_R3[j]; k++)
142  { // loop over stripeID
143  PHG4Hits.push_back(GetPHG4HitFromStripe(i, 3, j, k, electrons_per_stripe));
144  BotVertices.push_back(GetBotVerticesFromStripe(3, j, k));
145  TopVertices.push_back(GetTopVerticesFromStripe(3, j, k));
146  }
147  }
148  }
149 
150  // adjust G4Hits position and time
151  for (const auto& hit : PHG4Hits)
152  {
153  hit->set_t(0, m_centralMembraneDelay); //real hit delay
154  hit->set_t(1, m_centralMembraneDelay); //real hit delay.
155 
156  // TODO: check. There should be electrons on both sides on the central membrane
157  hit->set_z(0, 1.);
158  hit->set_z(1, 1.);
159  }
160 
162 }
163 
164 //_____________________________________________________________
166 {
167  // load g4hit container
168  auto g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
169  if (!g4hitcontainer)
170  {
171  std::cout << PHWHERE << "Could not locate g4 hit node " << hitnodename << std::endl;
173  }
174 
175  // copy all hits from G4hits vector into container
176  for (const auto& hit : PHG4Hits)
177  {
178  auto copy = new PHG4Hitv1(hit);
179  g4hitcontainer->AddHit(detId, copy);
180  }
181 
183 }
184 
185 //_____________________________________________________________
187 {
188  // same gas parameters as in PHG4TpcElectronDrift::SetDefaultParameters
189 
190  // Data on gasses @20 C and 760 Torr from the following source:
191  // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
192  // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
193  // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
194  static constexpr double Ne_dEdx = 1.56; // keV/cm
195  static constexpr double CF4_dEdx = 7.00; // keV/cm
196  static constexpr double Ne_NTotal = 43; // Number/cm
197  static constexpr double CF4_NTotal = 100; // Number/cm
198  static constexpr double Tpc_NTot = 0.5 * Ne_NTotal + 0.5 * CF4_NTotal;
199  static constexpr double Tpc_dEdx = 0.5 * Ne_dEdx + 0.5 * CF4_dEdx;
200  static constexpr double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
201 
202  // number of electrons per deposited GeV in TPC gas
203  set_default_double_param("electrons_per_gev", Tpc_ElectronsPerKeV * 1000000.);
204 
206  set_default_int_param("electrons_per_stripe", 300);
207 }
208 
209 //_____________________________________________________________
211  int nStripes, int nPads,
212  const std::array<double, nRadii>& R,
213  std::array<double, nRadii>& spacing,
214  double x1a[][nRadii], double y1a[][nRadii],
215  double x1b[][nRadii], double y1b[][nRadii],
216  double x2a[][nRadii], double y2a[][nRadii],
217  double x2b[][nRadii], double y2b[][nRadii],
218  double x3a[][nRadii], double y3a[][nRadii],
219  double x3b[][nRadii], double y3b[][nRadii],
220  double padfrac,
221  double str_width[][nRadii],
222  const std::array<double, nRadii>& widthmod,
223  std::array<int, nRadii>& nGoodStripes,
224  const std::array<int, nRadii>& keepUntil,
225  std::array<int, nRadii>& nStripesIn,
226  std::array<int, nRadii>& nStripesBefore)
227 {
228  const double phi_module = M_PI / 6.0; // angle span of a module
229  const int pr_mult = 3; // multiples of intrinsic resolution of pads
230  const int dw_mult = 8; // multiples of diffusion width
231  const double diffwidth = 0.6 * mm; // diffusion width
232  const double adjust = 0.015; //arbitrary angle to center the pattern in a petal
233 
234  double theta = 0.0;
235  //center coords
236  double cx[nStripes][nRadii], cy[nStripes][nRadii];
237  //corner coords
238  /* double tempX1a[nStripes][nRadii], tempY1a[nStripes][nRadii];
239  double tempX1b[nStripes][nRadii], tempY1b[nStripes][nRadii];
240  double tempX2a[nStripes][nRadii], tempY2a[nStripes][nRadii];
241  double tempX2b[nStripes][nRadii], tempY2b[nStripes][nRadii];
242  double rotatedX1a[nStripes][nRadii], rotatedY1a[nStripes][nRadii];
243  double rotatedX1b[nStripes][nRadii], rotatedY1b[nStripes][nRadii];
244  double rotatedX2a[nStripes][nRadii], rotatedY2a[nStripes][nRadii];
245  double rotatedX2b[nStripes][nRadii], rotatedY2b[nStripes][nRadii]; */
246 
247  //calculate spacing first:
248  for (int i = 0; i < nRadii; i++)
249  {
250  spacing[i] = 2.0 * ((dw_mult * diffwidth / R[i]) + (pr_mult * phi_module / nPads));
251  }
252 
253  //vertex calculation
254  for (int j = 0; j < nRadii; j++)
255  {
256  int i_out = 0;
257  for (int i = keepThisAndAfter[j]; i < keepUntil[j]; i++)
258  {
259  if (j % 2 == 0)
260  {
261  theta = i * spacing[j] + (spacing[j] / 2) - adjust;
262  cx[i_out][j] = R[j] * cos(theta);
263  cy[i_out][j] = R[j] * sin(theta);
264  }
265  else
266  {
267  theta = (i + 1) * spacing[j] - adjust;
268  cx[i_out][j] = R[j] * cos(theta);
269  cy[i_out][j] = R[j] * sin(theta);
270  }
271 
272  TVector3 corner[4];
273  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
274  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
275  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
276  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
277 
278  TVector3 rotatedcorner[4];
279  for (int i = 0; i < 4; i++)
280  {
281  rotatedcorner[i] = corner[i];
282  rotatedcorner[i].RotateZ(theta);
283  }
284 
285  x1a[i_out][j] = rotatedcorner[0].X() + cx[i_out][j];
286  x1b[i_out][j] = rotatedcorner[1].X() + cx[i_out][j];
287  x2a[i_out][j] = rotatedcorner[2].X() + cx[i_out][j];
288  x2b[i_out][j] = rotatedcorner[3].X() + cx[i_out][j];
289 
290  y1a[i_out][j] = rotatedcorner[0].Y() + cy[i_out][j];
291  y1b[i_out][j] = rotatedcorner[1].Y() + cy[i_out][j];
292  y2a[i_out][j] = rotatedcorner[2].Y() + cy[i_out][j];
293  y2b[i_out][j] = rotatedcorner[3].Y() + cy[i_out][j];
294 
295  /* x1a[i_out][j] = cx[i_out][j] - padfrac + arc_r;
296  y1a[i_out][j] = cy[i_out][j] - str_width/2;
297  x1b[i_out][j] = cx[i_out][j] + padfrac - arc_r;
298  y1b[i_out][j] = cy[i_out][j] - str_width/2;
299  x2a[i_out][j] = cx[i_out][j] - padfrac + arc_r;
300  y2a[i_out][j] = cy[i_out][j] + str_width/2;
301  x2b[i_out][j] = cx[i_out][j] + padfrac - arc_r;
302  y2b[i_out][j] = cy[i_out][j] + str_width/2;
303 
304  tempX1a[i_out][j] = x1a[i_out][j] - cx[i_out][j];
305  tempY1a[i_out][j] = y1a[i_out][j] - cy[i_out][j];
306  tempX1b[i_out][j] = x1b[i_out][j] - cx[i_out][j];
307  tempY1b[i_out][j] = y1b[i_out][j] - cy[i_out][j];
308  tempX2a[i_out][j] = x2a[i_out][j] - cx[i_out][j];
309  tempY2a[i_out][j] = y2a[i_out][j] - cy[i_out][j];
310  tempX2b[i_out][j] = x2b[i_out][j] - cx[i_out][j];
311  tempY2b[i_out][j] = y2b[i_out][j] - cy[i_out][j];
312 
313  rotatedX1a[i_out][j] = tempX1a[i_out][j]*cos(theta) - tempY1a[i_out][j]*sin(theta);
314  rotatedY1a[i_out][j] = tempX1a[i_out][j]*sin(theta) + tempY1a[i_out][j]*cos(theta);
315  rotatedX1b[i_out][j] = tempX1b[i_out][j]*cos(theta) - tempY1b[i_out][j]*sin(theta);
316  rotatedY1b[i_out][j] = tempX1b[i_out][j]*sin(theta) + tempY1b[i_out][j]*cos(theta);
317  rotatedX2a[i_out][j] = tempX2a[i_out][j]*cos(theta) - tempY2a[i_out][j]*sin(theta);
318  rotatedY2a[i_out][j] = tempX2a[i_out][j]*sin(theta) + tempY2a[i_out][j]*cos(theta);
319  rotatedX2b[i_out][j] = tempX2b[i_out][j]*cos(theta) - tempY2b[i_out][j]*sin(theta);
320  rotatedY2b[i_out][j] = tempX2b[i_out][j]*sin(theta) + tempY2b[i_out][j]*cos(theta);*/
321 
322  /* x1a[i_out][j] = rotatedX1a[i_out][j] + cx[i_out][j];
323  y1a[i_out][j] = rotatedY1a[i_out][j] + cy[i_out][j];
324  x1b[i_out][j] = rotatedX1b[i_out][j] + cx[i_out][j];
325  y1b[i_out][j] = rotatedY1b[i_out][j] + cy[i_out][j];
326  x2a[i_out][j] = rotatedX2a[i_out][j] + cx[i_out][j];
327  y2a[i_out][j] = rotatedY2a[i_out][j] + cy[i_out][j];
328  x2b[i_out][j] = rotatedX2b[i_out][j] + cx[i_out][j];
329  y2b[i_out][j] = rotatedY2b[i_out][j] + cy[i_out][j]; */
330 
331  x3a[i_out][j] = (x1a[i_out][j] + x2a[i_out][j]) / 2.0;
332  y3a[i_out][j] = (y1a[i_out][j] + y2a[i_out][j]) / 2.0;
333  x3b[i_out][j] = (x1b[i_out][j] + x2b[i_out][j]) / 2.0;
334  y3b[i_out][j] = (y1b[i_out][j] + y2b[i_out][j]) / 2.0;
335 
336  i_out++;
337 
338  nStripesBefore_R1_e[0] = 0;
339 
340  nStripesIn[j] = keepUntil[j] - keepThisAndAfter[j];
341 
342  nStripesBefore[j] = nStripesIn[j - 1] + nStripesBefore[j - 1];
343  nStripesBefore_R1_e[0] = 0;
344  }
345  nGoodStripes[j] = i_out;
346  }
347 }
348 
349 PHG4Hit* PHG4TpcCentralMembrane::GetBotVerticesFromStripe(int moduleID, int radiusID, int stripeID) const
350 {
351  PHG4Hit* botvert;
352  TVector3 dummyPos0, dummyPos1;
353 
354  //0 - left, 1 - right
355 
356  botvert = new PHG4Hitv1();
357  botvert->set_layer(-2);
358  if (moduleID == 0)
359  {
360  botvert->set_x(0, x1a_R1_e[stripeID][radiusID]);
361  botvert->set_y(0, y1a_R1_e[stripeID][radiusID]);
362  botvert->set_x(1, x1b_R1_e[stripeID][radiusID]);
363  botvert->set_y(1, y1b_R1_e[stripeID][radiusID]);
364  }
365  else if (moduleID == 1)
366  {
367  botvert->set_x(0, x1a_R1[stripeID][radiusID]);
368  botvert->set_y(0, y1a_R1[stripeID][radiusID]);
369  botvert->set_x(1, x1b_R1[stripeID][radiusID]);
370  botvert->set_y(1, y1b_R1[stripeID][radiusID]);
371  }
372  else if (moduleID == 2)
373  {
374  botvert->set_x(0, x1a_R2[stripeID][radiusID]);
375  botvert->set_y(0, y1a_R2[stripeID][radiusID]);
376  botvert->set_x(1, x1b_R2[stripeID][radiusID]);
377  botvert->set_y(1, y1b_R2[stripeID][radiusID]);
378  }
379  else if (moduleID == 3)
380  {
381  botvert->set_x(0, x1a_R3[stripeID][radiusID]);
382  botvert->set_y(0, y1a_R3[stripeID][radiusID]);
383  botvert->set_x(1, x1b_R3[stripeID][radiusID]);
384  botvert->set_y(1, y1b_R3[stripeID][radiusID]);
385  }
386  botvert->set_z(0, 0.0);
387 
388  return botvert;
389 }
390 
391 PHG4Hit* PHG4TpcCentralMembrane::GetTopVerticesFromStripe(int moduleID, int radiusID, int stripeID) const
392 {
393  PHG4Hit* topvert = nullptr;
394  TVector3 dummyPos0, dummyPos1;
395 
396  //0 - left, 1 - right
397 
398  topvert = new PHG4Hitv1();
399  topvert->set_layer(-3);
400  if (moduleID == 0)
401  {
402  topvert->set_x(0, x2a_R1_e[stripeID][radiusID]);
403  topvert->set_y(0, y2a_R1_e[stripeID][radiusID]);
404  topvert->set_x(1, x2b_R1_e[stripeID][radiusID]);
405  topvert->set_y(1, y2b_R1_e[stripeID][radiusID]);
406  }
407  else if (moduleID == 1)
408  {
409  topvert->set_x(0, x2a_R1[stripeID][radiusID]);
410  topvert->set_y(0, y2a_R1[stripeID][radiusID]);
411  topvert->set_x(1, x2b_R1[stripeID][radiusID]);
412  topvert->set_y(1, y2b_R1[stripeID][radiusID]);
413  }
414  else if (moduleID == 2)
415  {
416  topvert->set_x(0, x2a_R2[stripeID][radiusID]);
417  topvert->set_y(0, y2a_R2[stripeID][radiusID]);
418  topvert->set_x(1, x2b_R2[stripeID][radiusID]);
419  topvert->set_y(1, y2b_R2[stripeID][radiusID]);
420  }
421  else if (moduleID == 3)
422  {
423  topvert->set_x(0, x2a_R3[stripeID][radiusID]);
424  topvert->set_y(0, y2a_R3[stripeID][radiusID]);
425  topvert->set_x(1, x2b_R3[stripeID][radiusID]);
426  topvert->set_y(1, y2b_R3[stripeID][radiusID]);
427  }
428  topvert->set_z(0, 0.0);
429 
430  return topvert;
431 }
432 
434  const double x1a[][nRadii], const double x1b[][nRadii],
435  const double x2a[][nRadii], const double x2b[][nRadii],
436  const double y1a[][nRadii], const double y1b[][nRadii],
437  const double y2a[][nRadii], const double y2b[][nRadii],
438  const double x3a[][nRadii], const double y3a[][nRadii],
439  const double x3b[][nRadii], const double y3b[][nRadii],
440  double x, double y, const std::array<int, nRadii>& nGoodStripes) const
441 {
442  int c = 0;
443 
444  for (int j = 0; j < nRadii; j++)
445  {
446  for (int i = 0; i < nGoodStripes[j]; i++)
447  {
448  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])))
449  c = !c;
450  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])))
451  c = !c;
452  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])))
453  c = !c;
454  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])))
455  c = !c;
456 
457  //check inside arcs
458  if (c == 0)
459  {
460  if (((x - x3a[i][j]) * (x - x3a[i][j]) + (y - y3a[i][j]) * (y - y3a[i][j])) <= arc_r * arc_r)
461  {
462  c = !c;
463  }
464  else if (((x - x3b[i][j]) * (x - x3b[i][j]) + (y - y3b[i][j]) * (y - y3b[i][j])) <= arc_r * arc_r)
465  {
466  c = !c;
467  }
468  }
469  }
470  }
471  return c;
472 }
473 
474 int PHG4TpcCentralMembrane::getSearchResult(double xcheck, double ycheck) const
475 {
476  const double phi_petal = M_PI / 9.0; // angle span of one petal
477  const double end_R1_e = 312.0 * mm; // arbitrary radius between R1_e and R1
478  const double end_R1 = 408.0 * mm; // arbitrary radius between R1 and R2
479  const double end_R2 = 580.0 * mm; // arbitrary radius between R2 and R3
480 
481  double r, phi, phimod, xmod, ymod;
482 
483  r = sqrt(xcheck * xcheck + ycheck * ycheck);
484  phi = atan(ycheck / xcheck);
485  if ((xcheck < 0.0) && (ycheck > 0.0))
486  {
487  phi = phi + M_PI;
488  }
489  else if ((xcheck > 0.0) && (ycheck < 0.0))
490  {
491  phi = phi + 2.0 * M_PI;
492  }
493 
494  phimod = fmod(phi, phi_petal);
495  xmod = r * cos(phimod);
496  ymod = r * sin(phimod);
497 
498  int result = 0;
499 
500  if (r <= end_R1_e)
501  {
503  }
504  else if ((r > end_R1_e) && (r <= end_R1))
505  {
507  }
508  else if ((r > end_R1) && (r <= end_R2))
509  {
511  }
512  else if ((r > end_R2) && (r <= end_CM))
513  {
515  }
516 
517  return result;
518 }
519 
520 PHG4Hit* PHG4TpcCentralMembrane::GetPHG4HitFromStripe(int petalID, int moduleID, int radiusID, int stripeID, int nElectrons) const
521 { //this function generates a PHG4 hit using coordinates from a stripe
522  const double phi_petal = M_PI / 9.0; // angle span of one petal
523  PHG4Hit* hit;
524  TVector3 dummyPos0, dummyPos1;
525 
526  //could put in some sanity checks here but probably not necessary since this is only really used within the class
527  //petalID ranges 0-17, module ID 0-3, stripeID varies - nGoodStripes for each module
528  //radiusID ranges 0-7
529 
530  //from phg4tpcsteppingaction.cc
531  hit = new PHG4Hitv1();
532  hit->set_layer(-1); // dummy number
533  //here we set the entrance values in cm
534  if (moduleID == 0)
535  {
536  hit->set_x(0, x3a_R1_e[stripeID][radiusID] / cm);
537  hit->set_y(0, y3a_R1_e[stripeID][radiusID] / cm);
538  }
539  else if (moduleID == 1)
540  {
541  hit->set_x(0, x3a_R1[stripeID][radiusID] / cm);
542  hit->set_y(0, y3a_R1[stripeID][radiusID] / cm);
543  }
544  else if (moduleID == 2)
545  {
546  hit->set_x(0, x3a_R2[stripeID][radiusID] / cm);
547  hit->set_y(0, y3a_R2[stripeID][radiusID] / cm);
548  }
549  else if (moduleID == 3)
550  {
551  hit->set_x(0, x3a_R3[stripeID][radiusID] / cm);
552  hit->set_y(0, y3a_R3[stripeID][radiusID] / cm);
553  }
554  hit->set_z(0, 0.0 / cm);
555 
556  // check if you need to rotate coords to another petal
557  if (petalID > 0)
558  {
559  dummyPos0.SetXYZ(hit->get_x(0), hit->get_y(0), hit->get_z(0));
560  dummyPos0.RotateZ(petalID * phi_petal);
561  hit->set_x(0, dummyPos0.X());
562  hit->set_y(0, dummyPos0.Y());
563  }
564 
565  // TODO: use actual stripe direction for the momentum
566  hit->set_px(1, 500.0);
567  hit->set_py(1, 500.0);
568  hit->set_pz(1, 500.0);
569 
570  // time in ns
571  hit->set_t(0, 0);
572 
573  //set and save the track ID
574  hit->set_trkid(-1); // dummy number
575 
576  // here we just update the exit values, it will be overwritten
577  // for every step until we leave the volume or the particle
578  // ceases to exist
579  if (moduleID == 0)
580  {
581  hit->set_x(1, x3b_R1_e[stripeID][radiusID] / cm);
582  hit->set_y(1, y3b_R1_e[stripeID][radiusID] / cm);
583  }
584  else if (moduleID == 1)
585  {
586  hit->set_x(1, x3b_R1[stripeID][radiusID] / cm);
587  hit->set_y(1, y3b_R1[stripeID][radiusID] / cm);
588  }
589  else if (moduleID == 2)
590  {
591  hit->set_x(1, x3b_R2[stripeID][radiusID] / cm);
592  hit->set_y(1, y3b_R2[stripeID][radiusID] / cm);
593  }
594  else if (moduleID == 3)
595  {
596  hit->set_x(1, x3b_R3[stripeID][radiusID] / cm);
597  hit->set_y(1, y3b_R3[stripeID][radiusID] / cm);
598  }
599  hit->set_z(1, 0.0 / cm);
600 
601  // check if you need to rotate coords to another petal
602  if (petalID > 0)
603  {
604  dummyPos1.SetXYZ(hit->get_x(1), hit->get_y(1), hit->get_z(1));
605  dummyPos1.RotateZ(petalID * phi_petal);
606  hit->set_x(1, dummyPos1.X());
607  hit->set_y(1, dummyPos1.Y());
608  }
609 
610  // TODO: use actual stripe direction for the momentum
611  hit->set_px(1, 500.0);
612  hit->set_py(1, 500.0);
613  hit->set_pz(1, 500.0);
614 
615  hit->set_t(1, 0); // dummy number, nanosecond
616 
617  // calculate deposited energy corresponding to number of electrons per stripe
618  const double edep = nElectrons / electrons_per_gev;
619  hit->set_edep(edep);
620  hit->set_eion(edep);
621 
622  /*
623  if (hit->get_edep()){ //print out hits
624  double rin = sqrt(hit->get_x(0) * hit->get_x(0) + hit->get_y(0) * hit->get_y(0));
625  double rout = sqrt(hit->get_x(1) * hit->get_x(1) + hit->get_y(1) * hit->get_y(1));
626  std::cout << "Added Tpc g4hit with rin, rout = " << rin << " " << rout
627  << " g4hitid " << hit->get_hit_id() << std::endl;
628  std::cout << " xin " << hit->get_x(0)
629  << " yin " << hit->get_y(0)
630  << " zin " << hit->get_z(0)
631  << " rin " << rin
632  << std::endl;
633  std::cout << " xout " << hit->get_x(1)
634  << " yout " << hit->get_y(1)
635  << " zout " << hit->get_z(1)
636  << " rout " << rout
637  << std::endl;
638  std::cout << " xav " << (hit->get_x(1) + hit->get_x(0)) / 2.0
639  << " yav " << (hit->get_y(1) + hit->get_y(0)) / 2.0
640  << " zav " << (hit->get_z(1) + hit->get_z(0)) / 2.0
641  << " rav " << (rout + rin) / 2.0
642  << std::endl;
643  }
644  */
645 
646  return hit;
647 }
648 
649 int PHG4TpcCentralMembrane::getStripeID(double xcheck, double ycheck) const
650 {
651  //check if point came from stripe then see which stripe it is
652  //213 stripes in a petal, 18 petals, ntotstripes = 3834
653  int result, rID, petalID;
654  int phiID = 0;
655  int fullID = -1;
656  //double theta, spacing[nRadii], angle, m, dist;
657  double m, dist;
658  //const double adjust = 0.015; //arbitrary angle to center the pattern in a petal
659  const double phi_petal = M_PI / 9.0; // angle span of one petal
660 
661  double r, phi, phimod, xmod, ymod;
662 
663  // check if in a stripe
664  result = getSearchResult(xcheck, ycheck);
665 
666  // find which stripe
667  if (result == 1)
668  {
669  //std::cout << "on a stripe" << std::endl;
670  //convert coords to radius n angle
671  r = sqrt(xcheck * xcheck + ycheck * ycheck);
672  phi = atan(ycheck / xcheck);
673  if ((xcheck < 0.0) && (ycheck > 0.0))
674  {
675  phi = phi + M_PI;
676  }
677  else if ((xcheck > 0.0) && (ycheck < 0.0))
678  {
679  phi = phi + 2.0 * M_PI;
680  }
681  //get angle within first petal
682  phimod = fmod(phi, phi_petal);
683  xmod = r * cos(phimod);
684  ymod = r * sin(phimod);
685 
686  petalID = phi / phi_petal;
687 
688  for (int j = 0; j < nRadii; j++)
689  {
690  if (((R1_e[j] - padfrac_R1) < r) && (r < (R1_e[j] + padfrac_R1)))
691  { // check if radius is in stripe
692  rID = j;
693  std::cout << "rID: " << rID << std::endl;
694  //'angle' is to the center of a stripe
695  for (int i = 0; i < nGoodStripes_R1_e[j]; i++)
696  {
697  //if (j % 2 == 0){
698  //theta = i*spacing[j];
699  //angle = theta + (spacing[j]/2) - adjust;
700  // look at distance from center line of stripe
701  // if distance from x,y to center line < str_width
702  // calculate slope n then do dist
703 
704  m = (y3b_R1_e[i][j] - y3a_R1_e[i][j]) / (x3b_R1_e[i][j] - x3a_R1_e[i][j]);
705  /*std::cout << "y2: " << y3b_R1_e[i][j] << std::endl;
706  std::cout << "y1: " << y3a_R1_e[i][j] << std::endl;
707  std::cout << "x2: " << x3b_R1_e[i][j] << std::endl;
708  std::cout << "x1: " << x3a_R1_e[i][j] << std::endl;
709  std::cout << "xc: " << xcheck << std::endl;
710  std::cout << "yc: " << ycheck << std::endl;
711  std::cout << "m: " << m << std::endl; */
712  //std::cout << fabs((-m)*xcheck + ycheck) << std::endl;
713  dist = fabs((-m) * xmod + ymod) / sqrt(1 + m * m);
714  //std::cout << "dist:" << dist << std::endl;
715  if (dist < ((widthmod_R1_e[j] * str_width_R1_e[i][j]) / 2.0))
716  {
717  phiID = i;
718  //std::cout << "phiID: " << phiID << std::endl;
719  }
720  }
721 
722  std::cout << "nStripesBefore: " << nStripesBefore_R1_e[j] << std::endl;
723  fullID = petalID * nStripesPerPetal + nStripesBefore_R1_e[j] + phiID;
724  //std::cout << "fullID: " << fullID << std::endl;
725  }
726  else if (((R1[j] - padfrac_R1) < r) && (r < (R1[j] + padfrac_R1)))
727  {
728  rID = j + nRadii;
729  //std::cout << "R1" << std::endl;
730  for (int i = 0; i < nGoodStripes_R1[j]; i++)
731  {
732  // look at distance from center line of stripe
733  m = (y3b_R1[i][j] - y3a_R1[i][j]) / (x3b_R1[i][j] - x3a_R1[i][j]);
734  dist = fabs(m * xmod - ymod) / sqrt(1 + m * m);
735  if (dist < ((widthmod_R1[j] * str_width_R1[i][j]) / 2.0))
736  {
737  phiID = i;
738  }
739  }
740 
741  fullID = petalID * nStripesPerPetal + nStripesBefore_R1[j] + phiID;
742  }
743  else if (((R2[j] - padfrac_R2) < r) && (r < (R2[j] + padfrac_R2)))
744  {
745  rID = j + (2 * nRadii);
746  //std::cout << "R2" << std::endl;
747  for (int i = 0; i < nGoodStripes_R2[j]; i++)
748  {
749  // look at distance from center line of stripe
750  m = (y3b_R2[i][j] - y3a_R2[i][j]) / (x3b_R2[i][j] - x3a_R2[i][j]);
751  dist = fabs(m * xmod - ymod) / sqrt(1 + m * m);
752  if (dist < ((widthmod_R2[j] * str_width_R2[i][j]) / 2.0))
753  {
754  phiID = i;
755  }
756  }
757 
758  fullID = petalID * nStripesPerPetal + nStripesBefore_R2[j] + phiID;
759  }
760  else if (((R3[j] - padfrac_R3) < r) && (r < (R3[j] + padfrac_R3)))
761  {
762  rID = j + (3 * nRadii);
763  //std::cout << "R3" << std::endl;
764  for (int i = 0; i < nGoodStripes_R3[j]; i++)
765  {
766  // look at distance from center line of stripe
767  m = (y3b_R3[i][j] - y3a_R3[i][j]) / (x3b_R3[i][j] - x3a_R3[i][j]);
768  dist = fabs(m * xmod - ymod) / sqrt(1 + m * m);
769  if (dist < ((widthmod_R3[j] * str_width_R3[i][j]) / 2.0))
770  {
771  phiID = i;
772  }
773  }
774 
775  fullID = petalID * nStripesPerPetal + nStripesBefore_R3[j] + phiID;
776  }
777  }
778  }
779  else
780  {
781  fullID = -1;
782  }
783 
784  return fullID;
785 }