Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnnularFieldSim.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AnnularFieldSim.cc
1 #include "AnnularFieldSim.h"
2 
3 #include "AnalyticFieldModel.h"
4 #include "ChargeMapReader.h"
5 #include "MultiArray.h" //for TH3 alternative
6 #include "Rossegger.h"
7 
8 #include <TCanvas.h>
9 #include <TFile.h>
10 #include <TH1.h>
11 #include <TH2.h>
12 #include <TH3.h>
13 #include <TLatex.h>
14 #include <TPad.h>
15 #include <TString.h>
16 #include <TStyle.h>
17 #include <TTree.h>
18 #include <TVector3.h>
19 
20 #include <boost/format.hpp>
21 
22 #include <cassert> // for assert
23 #include <cmath>
24 #include <cstdio> // for printf, sprintf
25 #include <cstdlib> // for abs
26 #include <iostream>
27 
28 #define ALMOST_ZERO 0.00001
29 
30 AnnularFieldSim::AnnularFieldSim(float in_innerRadius, float in_outerRadius, float in_outerZ,
31  int r, int roi_r0, int roi_r1, int /*in_rLowSpacing*/, int /*in_rHighSize*/,
32  int phi, int roi_phi0, int roi_phi1, int /*in_phiLowSpacing*/, int /*in_phiHighSize*/,
33  int z, int roi_z0, int roi_z1, int /*in_zLowSpacing*/, int /*in_zHighSize*/,
34  float vdr, LookupCase in_lookupCase, ChargeCase in_chargeCase)
35 {
36  // check well-ordering:
37  if (roi_r0 >= r || roi_r1 > r || roi_r0 >= roi_r1)
38  {
39  assert(1 == 2);
40  }
41  if (roi_phi0 >= phi || roi_phi1 > phi || roi_phi0 >= roi_phi1)
42  {
43  printf("phi roi is out of range or spans the wrap-around. Please spare me that math.\n");
44  assert(1 == 2);
45  }
46  if (roi_z0 >= z || roi_z1 > z || roi_z0 >= roi_z1)
47  {
48  assert(1 == 2);
49  }
50 
51  hasTwin = false; // we never have a twin to start with.
52  green_shift = 0; // and so we don't shift our greens functions unless told to.
53 
54  printf("AnnularFieldSim::AnnularFieldSim with (%dx%dx%d) grid\n", r, phi, z);
55  printf("units m=%1.2E, cm=%1.2E, mm=%1.2E, um=%1.2E,\n", m, cm, mm, um);
56  printf("units s=%1.2E, us=%1.2E, ns=%1.2E,\n", s, us, ns);
57  printf("units C=%1.2E, nC=%1.2E, fC=%1.2E, \n", C, nC, fC);
58  printf("units Tesla=%1.2E, kGauss=%1.2E\n", Tesla, kGauss);
59 
60  // debug defaults:
61  //
64  debug_distortionScale.SetXYZ(0, 0, 0);
65  debug_npercent = 5;
66 
67  // internal states used for the debug flag
68  // goes with: if(debugFlag()) printf("%d: blah\n",__LINE__);
69 
70  // load constants of motion, dimensions, etc:
71  Enominal = 400 * V / cm; // v/cm
72  Bnominal = 1.4 * Tesla; // Tesla
73  vdrift = vdr * cm / s; // cm/s
74  langevin_T1 = 1.0;
75  langevin_T2 = 1.0;
76  omegatau_nominal = -10 * (Bnominal / kGauss) * (vdrift / (cm / us)) / (Enominal / (V / cm)); // to match to the familiar formula
77 
78  rmin = in_innerRadius * cm;
79  rmax = in_outerRadius * cm;
80  zmin = 0;
81  zmax = in_outerZ * cm;
82  if (zmin > zmax)
83  {
84  zmin = zmax;
85  zmax = 0; // for now, we continue to assume one end of the TPC is at z=0.
86  }
87  zero_vector.SetXYZ(0, 0, 0);
88 
89  // define the size of the volume:
90  dim.SetXYZ(1, 0, 0);
91  dim.SetPerp(rmax - rmin);
92  dim.SetPhi(0);
93  phispan = 2 * M_PI;
94  dim.SetZ(zmax - zmin);
95 
96  // set the green's functions model:
97  // UseFreeSpaceGreens();
98  // blah
99  green = nullptr;
100 
101  // load parameters of the whole-volume tiling
102  nr = r;
103  nphi = phi;
104  nz = z; // number of fundamental bins (f-bins) in each direction
105  printf("AnnularFieldSim::AnnularFieldSim set variables nr=%d, nphi=%d, nz=%d\n", nr, nphi, nz);
106  // and set the default truncation distance:
107  truncation_length = -1; // anything <1 and it won't truncate.
108 
109  // calculate the size of an f-bin:
110  // note that you have to set a non-zero value to start or perp won't update.
111  step.SetXYZ(1, 0, 0);
112  step.SetPerp(dim.Perp() / nr);
113  step.SetPhi(phispan / nphi);
114  step.SetZ(dim.Z() / nz);
115  printf("f-bin size: r=%f,phi=%f, z=%f, wanted %f,%f\n", step.Perp(), step.Phi(), (rmax - rmin) / nr, (phispan / nphi), (zmax - zmin) / nz);
116 
117  // create an array to store the charge in each f-bin
118  // q = new MultiArray<double>(nr, nphi, nz);
119  // q->SetAll(0);
120  q = new ChargeMapReader(nr, rmin, rmax, nphi, 0, phispan, nz, zmin, zmax);
121  sprintf(chargestring, "No spacecharge present.");
122 
123  // load parameters of our region of interest
124  rmin_roi = roi_r0;
125  phimin_roi = roi_phi0;
126  zmin_roi = roi_z0; // lower edge of our region of interest, measured in f-bins
127  rmax_roi = roi_r1;
128  phimax_roi = roi_phi1;
129  zmax_roi = roi_z1; // exlcuded upper edge of our region of interest, measured in f-bins
130  printf("AnnularFieldSim::AnnularFieldSim set roi variables rmin=%d phimin=%d zmin=%d rmax=%d phimax=%d zmax=%d\n",
132  // calculate the dimensions, in f-bins in our region of interest
136  printf("AnnularFieldSim::AnnularFieldSim calc'd roi variables nr=%d nphi=%d nz=%d\n", nr_roi, nphi_roi, nz_roi);
137 
138  // create an array to hold the SC-induced electric field in the roi with the specified dimensions
140  for (int i = 0; i < Efield->Length(); i++)
141  {
142  Efield->GetFlat(i)->SetXYZ(0, 0, 0);
143  }
144 
145  // and to hold the external electric fieldmap over the region of interest
147  for (int i = 0; i < Eexternal->Length(); i++)
148  {
149  Eexternal->GetFlat(i)->SetXYZ(0, 0, 0);
150  }
151 
152  // ditto the external magnetic fieldmap
154  for (int i = 0; i < Bfield->Length(); i++)
155  {
156  Bfield->GetFlat(i)->SetXYZ(0, 0, 0);
157  }
158 
159  // handle the lookup table construction:
160  lookupCase = in_lookupCase;
161  chargeCase = in_chargeCase;
162  if (chargeCase == ChargeCase::NoSpacecharge)
163  {
164  lookupCase = LookupCase::NoLookup; // don't build a lookup model if there's no charge. It just wastes time.
165  }
166 
167  if (lookupCase == Full3D)
168  {
169  printf("AnnularFieldSim::AnnularFieldSim building Epartial (full3D) with nr_roi=%d nphi_roi=%d nz_roi=%d =~%2.2fM TVector3 objects\n", nr_roi, nphi_roi, nz_roi,
170  nr_roi * nphi_roi * nz_roi * nr * nphi * nz / (1.0e6));
171 
173  for (int i = 0; i < Epartial->Length(); i++)
174  {
175  Epartial->GetFlat(i)->SetXYZ(0, 0, 0);
176  }
177  // and kill the arrays we shouldn't be using:
179  Epartial_highres->GetFlat(0)->SetXYZ(0, 0, 0);
180 
182  Epartial_lowres->GetFlat(0)->SetXYZ(0, 0, 0);
183 
185  Epartial_phislice->GetFlat(0)->SetXYZ(0, 0, 0);
186  q_lowres = new MultiArray<double>(1);
187  *(q_lowres->GetFlat(0)) = 0;
188  q_local = new MultiArray<double>(1);
189  *(q_local->GetFlat(0)) = 0;
190  }
191  else if (lookupCase == HybridRes)
192  {
193  printf("lookupCase==HybridRes\n");
194  // zero out the other two:
196  Epartial->GetFlat(0)->SetXYZ(0, 0, 0);
197 
199  Epartial_phislice->GetFlat(0)->SetXYZ(0, 0, 0);
200  }
201  else if (lookupCase == PhiSlice)
202  {
203  printf("lookupCase==PhiSlice\n");
204 
206  for (int i = 0; i < Epartial_phislice->Length(); i++)
207  {
208  Epartial_phislice->GetFlat(i)->SetXYZ(0, 0, 0);
209  }
210  // zero out the other two:
212  Epartial->GetFlat(0)->SetXYZ(0, 0, 0);
214  Epartial_highres->GetFlat(0)->SetXYZ(0, 0, 0);
215 
217  Epartial_lowres->GetFlat(0)->SetXYZ(0, 0, 0);
218 
219  q_lowres = new MultiArray<double>(1);
220  *(q_lowres->GetFlat(0)) = 0;
221  q_local = new MultiArray<double>(1);
222  *(q_local->GetFlat(0)) = 0;
223  }
224  else if (lookupCase == Analytic || lookupCase == NoLookup)
225  {
226  printf("lookupCase==Analytic (or NoLookup)\n");
227 
228  // zero them all out:
230  Epartial_phislice->GetFlat(0)->SetXYZ(0, 0, 0);
231 
233  Epartial->GetFlat(0)->SetXYZ(0, 0, 0);
234 
236  Epartial_highres->GetFlat(0)->SetXYZ(0, 0, 0);
237 
239  Epartial_lowres->GetFlat(0)->SetXYZ(0, 0, 0);
240 
241  q_lowres = new MultiArray<double>(1);
242  *(q_lowres->GetFlat(0)) = 0;
243  q_local = new MultiArray<double>(1);
244  *(q_local->GetFlat(0)) = 0;
245  }
246  else
247  {
248  printf("Ran into wrong lookupCase logic in constructor.\n");
249  assert(1 == 2);
250  }
251 
252  return;
253 }
254 AnnularFieldSim::AnnularFieldSim(float in_innerRadius, float in_outerRadius, float in_outerZ,
255  int r, int roi_r0, int roi_r1, int in_rLowSpacing, int in_rHighSize,
256  int phi, int roi_phi0, int roi_phi1, int in_phiLowSpacing, int in_phiHighSize,
257  int z, int roi_z0, int roi_z1, int in_zLowSpacing, int in_zHighSize,
258  float vdr, LookupCase in_lookupCase)
259  : AnnularFieldSim(in_innerRadius, in_outerRadius, in_outerZ,
260  r, roi_r0, roi_r1, in_rLowSpacing, in_rHighSize,
261  phi, roi_phi0, roi_phi1, in_phiLowSpacing, in_phiHighSize,
262  z, roi_z0, roi_z1, in_zLowSpacing, in_zHighSize,
263  vdr, in_lookupCase, ChargeCase::FromFile)
264 {
265  printf("AnnularFieldSim::OldConstructor building AnnularFieldSim with ChargeCase::FromFile\n");
266  return;
267 }
268 AnnularFieldSim::AnnularFieldSim(float in_innerRadius, float in_outerRadius, float in_outerZ,
269  int r, int roi_r0, int roi_r1,
270  int phi, int roi_phi0, int roi_phi1,
271  int z, int roi_z0, int roi_z1,
272  float vdr, LookupCase in_lookupCase)
273  : AnnularFieldSim(in_innerRadius, in_outerRadius, in_outerZ,
274  r, roi_r0, roi_r1, 1, 3,
275  phi, roi_phi0, roi_phi1, 1, 3,
276  z, roi_z0, roi_z1, 1, 3,
277  vdr, in_lookupCase, ChargeCase::FromFile)
278 {
279  // just passing through for the old version again
280  // creates a region with high-res size of 3 (enough to definitely cover the eight local l-bins) and low-res spacing of 1, which ought to match the behavior (with a little more overhead) from when there was no highres-lowres distinction
281  printf("AnnularFieldSim::OldConstructor building AnnularFieldSim with local_size=1 in all dimensions, lowres_spacing=1 in all dimensions\n");
282  return;
283 }
284 AnnularFieldSim::AnnularFieldSim(float rin, float rout, float dz, int r, int phi, int z, float vdr)
285  : AnnularFieldSim(rin, rout, dz,
286  r, 0, r,
287  phi, 0, phi,
288  z, 0, z,
289  vdr, LookupCase::PhiSlice)
290 {
291  // just a pass-through to go from the old style to the more detailed version.
292  printf("AnnularFieldSim::OldConstructor building AnnularFieldSim with roi=full in all dimensions\n");
293  return;
294 }
295 
296 TVector3 AnnularFieldSim::calc_unit_field(TVector3 at, TVector3 from)
297 {
298  // if(debugFlag()) printf("%d: AnnularFieldSim::calc_unit_field(at=(r=%f,phi=%f,z=%f))\n",__LINE__,at.Perp(),at.Phi(),at.Z());
299  int r_position;
300  if (GetRindexAndCheckBounds(at.Perp(), &r_position) != InBounds)
301  {
302  printf("something's asking for 'at' with r=%f, which is index=%d\n", at.Perp(), r_position);
303  assert(1 == 2);
304  }
305  // note this is the field due to a fixed point charge in free space.
306  // if doing cylindrical calcs with different boundary conditions, this needs to change.
307 
308  // this could check roi bounds before returning, if things start acting funny.
309 
310  TVector3 field(0, 0, 0);
311 
312  // if we're more than truncation_length away, don't bother? rcc food for thought. right now _length is in bins, so tricky.
313 
314  // if the green's function class is not present use free space:
315  if (green == nullptr)
316  {
317  TVector3 delr = at - from;
318  field = delr; // to set the direction.
319  if (delr.Mag() < ALMOST_ZERO * ALMOST_ZERO)
320  { // note that this has blurred units -- it should scale with all three dimensions of stepsize. For lots of phi bins, especially, this might start to read as small before it's really small.
321  // do nothing. the vector is already zero, which will be our approximation.
322  // field.SetMag(0);//no contribution if we're in the same cell. -- but root warns if trying to resize something of magnitude zero.
323  }
324  else
325  {
326  field.SetMag(k_perm * 1 / (delr * delr)); // scalar product on the bottom. unitful, since we defined k_perm and delr with their correct units. (native units V=1 C=1 cm=1)
327  }
328  // printf("calc_unit_field at (%2.2f,%2.2f,%2.2f) from (%2.2f,%2.2f,%2.2f). Mag=%2.4fe-9\n",at.x(),at.Y(),at.Z(),from.X(),from.Y(),from.Z(),field.Mag()*1e9);
329  }
330  else
331  {
332  double atphi = FilterPhiPos(at.Phi());
333  double fromphi = FilterPhiPos(from.Phi());
334  double delphi = abs(at.DeltaPhi(from));
335  // to allow us to use the same greens function set for both sides of the tpc, we shift into the valid greens region if needed:
336  at.SetZ(at.Z() + green_shift);
337  from.SetZ(from.Z() + green_shift);
338  double Er = green->Er(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
339  // RCC manually disabled phi component of green -- actually, a correction to disallow trying to compute phi terms when at the same phi:
340  double Ephi = 0;
341  if (delphi > ALMOST_ZERO)
342  {
343  Ephi = green->Ephi(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
344  }
345  double Ez = green->Ez(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
346  field.SetXYZ(-Er, -Ephi, -Ez); // now these are the correct components if our test point is at y=0 (hence phi=0);
347  field = field * epsinv; // scale field strength, since the greens functions as of Apr 1 2020 do not build-in this factor.
348  field.RotateZ(at.Phi()); // rotate to the coordinates of our 'at' point, which is a small rotation for the phislice case.
349  // but this does mean we need to be careful! If we are rotating so that this is pointing not in the R direction, then with azimuthally symmetric charge we will have some cancellation we don't want, maybe? Make sure this matches how we rotate to sum_field!
350  }
351  return field;
352 }
353 
355 {
356  // this primarily takes the region [-pi,0] and maps it to [pi,2pi] by adding 2pi to it.
357  // if math has pushed us past 2pi, it also subtracts to try to get us in range.
358  double p = phi;
359  if (p >= 2 * M_PI) // phispan)
360  {
361  p -= 2 * M_PI;
362  }
363  if (p < 0)
364  {
365  p += 2 * M_PI;
366  }
367  if (p >= 2 * M_PI || p < 0)
368  {
369  printf("AnnularFieldSim::FilterPhiPos asked to filter %f, which is more than range=%f out of bounds. Check what called this.\n", phi, 2 * M_PI);
370  assert(1 == 2);
371  }
372  return p;
373 }
374 int AnnularFieldSim::FilterPhiIndex(int phi, int range = -1)
375 {
376  if (range < 0)
377  {
378  range = nphi; // default input is range=-1. in that case, use the intrinsic resolution of the q grid.
379  }
380 
381  // shifts phi up or down by nphi (=2pi in phi index space), and complains if this doesn't put it in range.
382  int p = phi;
383  if (p >= range)
384  {
385  p -= range;
386  }
387  if (p < 0)
388  {
389  p += range;
390  }
391  if (p >= range || p < 0)
392  {
393  printf("AnnularFieldSim::FilterPhiIndex asked to filter %d, which is more than range=%d out of bounds. Check what called this.\n", phi, range);
394  assert(1 == 2);
395  }
396  return p;
397 }
398 
400 {
401  float r0f = (pos - rmin) / step.Perp(); // the position in r, in units of step, starting from the low edge of the 0th bin.
402  int r0 = floor(r0f);
403  return r0;
404 }
405 
407 {
408  float p0f = (pos) / step.Phi(); // the position in phi, in units of step, starting from the low edge of the 0th bin.
409  int phitemp = floor(p0f);
410  int p0 = FilterPhiIndex(phitemp);
411  return p0;
412 }
413 
415 {
416  float z0f = (pos - zmin) / step.Z(); // the position in z, in units of step, starting from the low edge of the 0th bin.
417  int z0 = floor(z0f);
418  return z0;
419 }
420 
422 {
423  // if(debugFlag()) printf("%d: AnnularFieldSim::GetRindexAndCheckBounds(r=%f)\n",__LINE__,pos);
424 
425  float r0f = (pos - rmin) / step.Perp(); // the position in r, in units of step, starting from the low edge of the 0th bin.
426  int r0 = floor(r0f);
427  *r = r0;
428 
429  int r0lowered_slightly = floor(r0f - ALMOST_ZERO);
430  int r0raised_slightly = floor(r0f + ALMOST_ZERO);
431  if (r0lowered_slightly >= rmax_roi || r0raised_slightly < rmin_roi)
432  {
433  return OutOfBounds;
434  }
435 
436  // now if we are out of bounds, it is because we are on an edge, within range of ALMOST_ZERO of being in fair territory.
437  if (r0 >= rmax_roi)
438  {
439  return OnHighEdge;
440  }
441  if (r0 < rmin_roi)
442  {
443  return OnLowEdge;
444  }
445  // if we're still here, we're in bounds.
446  return InBounds;
447 }
449 {
450  // if(debugFlag()) printf("%d: AnnularFieldSim::GetPhiIndexAndCheckBounds(phi=%f)\n\n",__LINE__,pos);
451  float p0f = (pos) / step.Phi(); // the position in phi, in units of step, starting from the low edge of the 0th bin.
452  int phitemp = floor(p0f);
453  int p0 = FilterPhiIndex(phitemp);
454  *phi = p0;
455 
456  phitemp = floor(p0f - ALMOST_ZERO);
457  int p0lowered_slightly = FilterPhiIndex(phitemp);
458  phitemp = floor(p0f + ALMOST_ZERO);
459  int p0raised_slightly = FilterPhiIndex(phitemp);
460  // annoying detail: if we are at index 0, we might go above pmax by going down.
461  // and if we are at nphi-1, we might go below pmin by going up.
462  // if we are not at p0=0 or nphi-1, the slight wiggles won't move us.
463  // if we are at p0=0, we are not at or above the max, and lowering slightly won't change that,
464  // and is we are at p0=nphi-1, we are not below the min, and raising slightly won't change that
465  if ((p0lowered_slightly >= phimax_roi && p0 != 0) || (p0raised_slightly < phimin_roi && p0 != (nphi - 1)))
466  {
467  return OutOfBounds;
468  }
469  // now if we are out of bounds, it is because we are on an edge, within range of ALMOST_ZERO of being in fair territory.
470  if (p0 >= phimax_roi)
471  {
472  return OnHighEdge;
473  }
474  if (p0 < phimin_roi)
475  {
476  return OnLowEdge;
477  }
478  // if we're still here, we're in bounds.
479  return InBounds;
480 }
481 
483 {
484  // if(debugFlag()) printf("%d: AnnularFieldSim::GetZindexAndCheckBounds(z=%f)\n\n",__LINE__,pos);
485  float z0f = (pos - zmin) / step.Z(); // the position in z, in units of step, starting from the low edge of the 0th bin.
486  int z0 = floor(z0f);
487  *z = z0;
488 
489  int z0lowered_slightly = floor(z0f - ALMOST_ZERO);
490  int z0raised_slightly = floor(z0f + ALMOST_ZERO);
491 
492  if (z0lowered_slightly >= zmax_roi || z0raised_slightly < zmin_roi)
493  {
494  return OutOfBounds;
495  }
496  // now if we are out of bounds, it is because we are on an edge, within range of ALMOST_ZERO of being in fair territory.
497  if (z0 >= zmax_roi)
498  {
499  return OnHighEdge;
500  }
501  if (z0 < zmin_roi)
502  {
503  return OnLowEdge;
504  }
505  // if we're still here, we're in bounds.
506  return InBounds;
507 }
508 
510 {
511  // integrates E dz, from the starting point to the selected z position. The path is assumed to be along z for each step, with adjustments to x and y accumulated after each step.
512  // if(debugFlag()) printf("%d: AnnularFieldSim::fieldIntegral(x=%f,y=%f, z=%f) to z=%f\n\n",__LINE__,start.X(),start.Y(),start.Z(),zdest);
513  // printf("AnnularFieldSim::analyticFieldIntegral calculating from (%f,%f,%f) (rphiz)=(%f,%f,%f) to z=%f.\n",start.X(),start.Y(),start.Z(),start.Perp(),start.Phi(),start.Z(),zdest);
514 
515  // coordinates are assumed to be in native units (cm=1);
516 
517  int r, phi;
518  bool rOkay = (GetRindexAndCheckBounds(start.Perp(), &r) == InBounds);
519  bool phiOkay = (GetPhiIndexAndCheckBounds(FilterPhiPos(start.Phi()), &phi) == InBounds);
520 
521  // bool isE=(field==Efield);
522  // bool isB=(field==Bfield);
523 
524  // printf("anaFieldInt: isE=%d, isB=%d, start=(%f,%f,%f)\n",isE,isB,start.X(),start.Y(),start.Z());
525 
526  if (!rOkay || !phiOkay)
527  {
528  printf("AnnularFieldSim::analyticFieldIntegral asked to integrate along (r=%f,phi=%f), index=(%d,%d), which is out of bounds. returning starting position.\n", start.Perp(), start.Phi(), r, phi);
529  return start;
530  }
531 
532  int dir = (start.Z() > zdest ? -1 : 1); //+1 if going to larger z, -1 if going to smaller; if they're the same, the sense doesn't matter.
533 
534  int zi, zf;
535  double startz, endz;
536  BoundsCase startBound, endBound;
537 
538  // make sure 'zi' is always the smaller of the two numbers, for handling the partial-steps.
539  if (dir > 0)
540  {
541  startBound = GetZindexAndCheckBounds(start.Z(), &zi); // highest cell with lower bound less than lower bound of integral
542  endBound = GetZindexAndCheckBounds(zdest, &zf); // highest cell with lower bound less than upper lower bound of integral
543  startz = start.Z();
544  endz = zdest;
545  }
546  else
547  {
548  endBound = GetZindexAndCheckBounds(start.Z(), &zf); // highest cell with lower bound less than lower bound of integral
549  startBound = GetZindexAndCheckBounds(zdest, &zi); // highest cell with lower bound less than upper lower bound of integral
550  startz = zdest;
551  endz = start.Z();
552  }
553  bool startOkay = (startBound == InBounds);
554  bool endOkay = (endBound == InBounds || endBound == OnHighEdge); // if we just barely touch out-of-bounds on the high end, we can skip that piece of the integral
555 
556  if (!startOkay || !endOkay)
557  {
558  printf("AnnularFieldSim::analyticFieldIntegral asked to integrate from z=%f to %f, index=%d to %d), which is out of bounds. returning starting position.\n", startz, endz, zi, zf);
559  return start;
560  }
561 
562  start.SetZ(startz);
563  TVector3 integral;
564  if (field == Efield)
565  {
566  integral = aliceModel->Eint(endz, start) + Eexternal->Get(r - rmin_roi, phi - phimin_roi, zi - zmin_roi) * (endz - startz);
567  return dir * integral;
568  }
569  else if (field == Bfield)
570  {
571  return interpolatedFieldIntegral(zdest, start, Bfield);
572  }
573  return integral;
574 }
575 
576 TVector3 AnnularFieldSim::fieldIntegral(float zdest, const TVector3 &start, MultiArray<TVector3> *field)
577 {
578  // integrates E dz, from the starting point to the selected z position. The path is assumed to be along z for each step, with adjustments to x and y accumulated after each step.
579  // if(debugFlag()) printf("%d: AnnularFieldSim::fieldIntegral(x=%f,y=%f, z=%f) to z=%f\n\n",__LINE__,start.X(),start.Y(),start.Z(),zdest);
580 
581  int r, phi;
582  bool rOkay = (GetRindexAndCheckBounds(start.Perp(), &r) == InBounds);
583  bool phiOkay = (GetPhiIndexAndCheckBounds(FilterPhiPos(start.Phi()), &phi) == InBounds);
584 
585  if (!rOkay || !phiOkay)
586  {
587  printf("AnnularFieldSim::fieldIntegral asked to integrate along (r=%f,phi=%f), index=(%d,%d), which is out of bounds. returning starting position.\n", start.Perp(), start.Phi(), r, phi);
588  return start;
589  }
590 
591  int dir = (start.Z() > zdest ? -1 : 1); //+1 if going to larger z, -1 if going to smaller; if they're the same, the sense doesn't matter.
592 
593  int zi, zf;
594  double startz, endz;
595  BoundsCase startBound, endBound;
596 
597  // make sure 'zi' is always the smaller of the two numbers, for handling the partial-steps.
598  if (dir > 0)
599  {
600  startBound = GetZindexAndCheckBounds(start.Z(), &zi); // highest cell with lower bound less than lower bound of integral
601  endBound = GetZindexAndCheckBounds(zdest, &zf); // highest cell with lower bound less than upper lower bound of integral
602  startz = start.Z();
603  endz = zdest;
604  }
605  else
606  {
607  endBound = GetZindexAndCheckBounds(start.Z(), &zf); // highest cell with lower bound less than lower bound of integral
608  startBound = GetZindexAndCheckBounds(zdest, &zi); // highest cell with lower bound less than upper lower bound of integral
609  startz = zdest;
610  endz = start.Z();
611  }
612  bool startOkay = (startBound == InBounds);
613  bool endOkay = (endBound == InBounds || endBound == OnHighEdge); // if we just barely touch out-of-bounds on the high end, we can skip that piece of the integral
614 
615  if (!startOkay || !endOkay)
616  {
617  printf("AnnularFieldSim::fieldIntegral asked to integrate from z=%f to %f, index=%d to %d), which is out of bounds. returning starting position.\n", startz, endz, zi, zf);
618  return start;
619  }
620 
621  TVector3 fieldInt(0, 0, 0);
622  // printf("AnnularFieldSim::fieldIntegral requesting (%d,%d,%d)-(%d,%d,%d) (inclusive) cells\n",r,phi,zi,r,phi,zf-1);
623  TVector3 tf;
624  for (int i = zi; i < zf; i++)
625  { // count the whole cell of the lower end, and skip the whole cell of the high end.
626  tf = field->Get(r - rmin_roi, phi - phimin_roi, i - zmin_roi);
627  // printf("fieldAt (%d,%d,%d)=(%f,%f,%f) step=%f\n",r,phi,i,tf.X(),tf.Y(),tf.Z(),step.Z());
628  fieldInt += tf * step.Z();
629  }
630 
631  // since bins contain their lower bound, but not their upper, I can safely remove the unused portion of the lower cell:
632  fieldInt -= field->Get(r - rmin_roi, phi - phimin_roi, zi - zmin_roi) * (startz - zi * step.Z()); // remove the part of the low end cell we didn't travel through
633 
634  // but only need to add the used portion of the upper cell if we go past the edge of it meaningfully:
635  if (endz / step.Z() - zf > ALMOST_ZERO)
636  {
637  // printf("endz/step.Z()=%f, zf=%f\n",endz/step.Z(),zf*1.0);
638  // if our final step is actually in the next step.
639  fieldInt += field->Get(r - rmin_roi, phi - phimin_roi, zf - zmin_roi) * (endz - zf * step.Z()); // add the part of the high end cell we did travel through
640  }
641 
642  return dir * fieldInt;
643 }
644 
645 TVector3 AnnularFieldSim::GetCellCenter(int r, int phi, int z)
646 {
647  // returns the midpoint of the cell (halfway between each edge, not weighted center)
648 
649  TVector3 c(1, 0, 0);
650  c.SetPerp((r + 0.5) * step.Perp() + rmin);
651  c.SetPhi((phi + 0.5) * step.Phi());
652  c.SetZ((z + 0.5) * step.Z() + zmin);
653 
654  return c;
655 }
656 
657 TVector3 AnnularFieldSim::GetRoiCellCenter(int r, int phi, int z)
658 {
659  // returns the midpoint of the cell relative to the roi (halfway between each edge, not weighted center)
660 
661  // return zero if it's out of bounds:
662  if (r > nr_roi || r < 0 || phi > nphi_roi || phi < 0 || z > nz_roi || z < 0)
663  {
664  return zero_vector;
665  }
666 
667  TVector3 c(1, 0, 0);
668  c.SetPerp((r + rmin_roi + 0.5) * step.Perp() + rmin);
669  c.SetPhi((phi + phimin_roi + 0.5) * step.Phi());
670  c.SetZ((z + zmin_roi + 0.5) * step.Z() + zmin);
671 
672  return c;
673 }
674 
675 TVector3 AnnularFieldSim::GetGroupCellCenter(int r0, int r1, int phi0, int phi1, int z0, int z1)
676 {
677  // returns the midpoint of the cell (halfway between each edge, not weighted center)
678  float ravg = (r0 + r1) / 2.0 + 0.5;
679  if (phi0 > phi1)
680  {
681  phi1 += nphi;
682  }
683  if (phi0 > phi1)
684  {
685  printf("phi1(%d)<=phi0(%d) even after boosting phi1. check what called this!\n", phi1, phi0);
686  assert(1 == 2);
687  }
688  float phiavg = (r0 + r1) / 2.0 + 0.5;
689  if (phiavg >= nphi)
690  {
691  phiavg -= nphi;
692  }
693 
694  float zavg = (z0 + z1) / 2.0 + 0.5;
695 
696  TVector3 c(1, 0, 0);
697  c.SetPerp((ravg) *step.Perp() + rmin);
698  c.SetPhi((phiavg) *step.Phi());
699  c.SetZ((zavg) *step.Z() + zmin);
700 
701  return c;
702 }
703 
705 {
706  // returns the weighted center of the cell by volume.
707  // todo: this is vaguely hefty, and might be worth storing the result of, if speed is needed
708  TVector3 c(1, 0, 0);
709 
710  float rin = r * step.Perp() + rmin;
711  float rout = rin + step.Perp();
712 
713  float rMid = (4 * sin(step.Phi() / 2) * (pow(rout, 3) - pow(rin, 3)) / (3 * step.Phi() * (pow(rout, 2) - pow(rin, 2))));
714  c.SetPerp(rMid);
715  c.SetPhi((phi + 0.5) * step.Phi());
716  c.SetZ((z + 0.5) * step.Z());
717 
718  return c;
719 }
720 
722 {
723  // printf("AnnularFieldSim::interpolatedFieldIntegral(x=%f,y=%f, z=%f)\n",start.X(),start.Y(),start.Z());
724 
725  float r0 = (start.Perp() - rmin) / step.Perp() - 0.5; // the position in r, in units of step, starting from the center of the 0th bin.
726  int r0i = floor(r0); // the integer portion of the position. -- what center is below our position?
727  float r0d = r0 - r0i; // the decimal portion of the position. -- how far past center are we?
728  int ri[4]; // the r position of the four cell centers to consider
729  ri[0] = ri[1] = r0i;
730  ri[2] = ri[3] = r0i + 1;
731  float rw[4]; // the weight of that cell, linear in distance from the center of it
732  rw[0] = rw[1] = 1 - r0d; // 1 when we're on it, 0 when we're at the other one.
733  rw[2] = rw[3] = r0d; // 1 when we're on it, 0 when we're at the other one
734 
735  bool skip[] = {false, false, false, false};
736  if (ri[0] < rmin_roi)
737  {
738  skip[0] = true; // don't bother handling 0 and 1 in the coordinates.
739  skip[1] = true;
740  rw[2] = rw[3] = 1; // and weight like we're dead-center on the outer cells.
741  }
742  if (ri[2] >= rmax_roi)
743  {
744  skip[2] = true; // don't bother handling 2 and 3 in the coordinates.
745  skip[3] = true;
746  rw[0] = rw[1] = 1; // and weight like we're dead-center on the inner cells.
747  }
748 
749  // now repeat that structure for phi:
750  float p0 = (start.Phi()) / step.Phi() - 0.5; // the position in phi, in units of step, starting from the center of the 0th bin.
751  int p0i = floor(p0); // the integer portion of the position. -- what center is below our position?
752  float p0d = p0 - p0i; // the decimal portion of the position. -- how far past center are we?
753  int pi[4]; // the phi position of the four cell centers to consider
754  pi[0] = pi[2] = FilterPhiIndex(p0i);
755  pi[1] = pi[3] = FilterPhiIndex(p0i + 1);
756  float pw[4]; // the weight of that cell, linear in distance from the center of it
757  pw[0] = pw[2] = 1 - p0d; // 1 when we're on it, 0 when we're at the other one.
758  pw[1] = pw[3] = p0d; // 1 when we're on it, 0 when we're at the other one
759 
760  // to handle wrap-around
761  if (pi[0] < phimin_roi || pi[0] >= phimax_roi)
762  {
763  skip[0] = true; // don't bother handling 0 and 1 in the coordinates.
764  skip[2] = true;
765  pw[1] = pw[3] = 1; // and weight like we're dead-center on the outer cells.
766  }
767  if (pi[1] >= phimax_roi || pi[1] < phimin_roi)
768  {
769  skip[1] = true; // don't bother handling 2 and 3 in the coordinates.
770  skip[3] = true;
771  pw[0] = pw[2] = 1; // and weight like we're dead-center on the inner cells.
772  }
773 
774  int dir = (start.Z() < zdest ? 1 : -1); //+1 if going to larger z, -1 if going to smaller; if they're the same, the sense doesn't matter.
775 
776  int zi, zf;
777  double startz, endz;
778  BoundsCase startBound, endBound;
779 
780  // make sure 'zi' is always the smaller of the two numbers, for handling the partial-steps.
781  if (dir > 0)
782  {
783  startBound = GetZindexAndCheckBounds(start.Z(), &zi); // highest cell with lower bound less than lower bound of integral
784  endBound = GetZindexAndCheckBounds(zdest, &zf); // highest cell with lower bound less than upper bound of integral
785  startz = start.Z();
786  endz = zdest;
787  }
788  else
789  {
790  endBound = GetZindexAndCheckBounds(start.Z(), &zf); // highest cell with lower bound less than lower bound of integral
791  startBound = GetZindexAndCheckBounds(zdest, &zi); // highest cell with lower bound less than upper bound of integral
792  startz = zdest;
793  endz = start.Z();
794  }
795  bool startOkay = (startBound == InBounds || startBound == OnLowEdge); // maybe todo: add handling for being just below the low edge.
796  bool endOkay = (endBound == InBounds || endBound == OnHighEdge); // if we just barely touch out-of-bounds on the high end, we can skip that piece of the integral
797 
798  if (!startOkay || !endOkay)
799  {
800  printf("AnnularFieldSim::InterpolatedFieldIntegral asked to integrate from z=%f to %f, index=%d to %d), which is out of bounds. returning zero\n", startz, endz, zi, zf);
801  return zero_vector;
802  }
803 
804  if (startBound == OnLowEdge)
805  {
806  // we were just below the low edge, so we will be asked to sample a bin in z we're not actually using
807  zi++; // avoid it. We weren't integrating anything in it anyway.
808  }
809 
810  if (false)
811  { // debugging options
812  bool isE = (field == Efield);
813  bool isB = (field == Bfield);
814  printf("interpFieldInt: isE=%d, isB=%d, start=(%2.4f,%2.4f,%2.4f) dest=%2.4f\n", isE, isB, start.X(), start.Y(), start.Z(), zdest);
815  printf("interpolating fieldInt for %s at r=%f,phi=%f\n", (field == Efield) ? "Efield" : "Bfield", r0, p0);
816  printf("direction=%d, startz=%2.4f, endz=%2.4f, zi=%d, zf=%d\n", dir, startz, endz, zi, zf);
817  for (int i = 0; i < 4; i++)
818  {
819  printf("skip[%d]=%d\tw[i]=%2.2f, (pw=%2.2f, rw=%2.2f)\n", i, skip[i], rw[i] * pw[i], pw[i], rw[i]);
820  }
821  }
822 
823  TVector3 fieldInt(0, 0, 0), partialInt; // where we'll store integrals as we generate them.
824 
825  for (int i = 0; i < 4; i++)
826  {
827  if (skip[i])
828  {
829  // printf("skipping element r=%d,phi=%d\n",ri[i],pi[i]);
830  continue; // we invalidated this one for some reason.
831  }
832  partialInt.SetXYZ(0, 0, 0);
833  for (int j = zi; j < zf; j++)
834  { // count the whole cell of the lower end, and skip the whole cell of the high end.
835 
836  partialInt += field->Get(ri[i] - rmin_roi, pi[i] - phimin_roi, j - zmin_roi) * step.Z();
837  }
838  if (startBound != OnLowEdge)
839  {
840  partialInt -= field->Get(ri[i] - rmin_roi, pi[i] - phimin_roi, zi - zmin_roi) * (startz - (zi * step.Z() + zmin)); // remove the part of the low end cell we didn't travel through
841  // printf("removing low end of cell we didn't travel through (zi-zmin_roi=%d, length=%f)\n",zi-zmin_roi, startz-(zi*step.Z()+zmin));
842  }
843  if ((endz - zmin) / step.Z() - zf > ALMOST_ZERO)
844  {
845  partialInt += field->Get(ri[i] - rmin_roi, pi[i] - phimin_roi, zf - zmin_roi) * (endz - (zf * step.Z() + zmin)); // add the part of the high end cell we did travel through
846  // printf("adding low end of cell we did travel through (zf-zmin_roi=%d, length=%f)\n",zf-zmin_roi, endz-(zf*step.Z()+zmin));
847  }
848  // printf("element r=%d,phi=%d, w=%f partialInt=(%2.2E,%2.2E,%2.2E)\n",ri[i],pi[i],rw[i]*pw[i],partialInt.X(),partialInt.Y(),partialInt.Z());
849 
850  fieldInt += rw[i] * pw[i] * partialInt;
851  }
852 
853  return dir * fieldInt;
854 }
855 
857 {
858  // scalefactor should be chosen so Rho(pos) returns C/cm^3
859 
860  // sphenix:
861  // double ifc_radius=20;
862  // double ofc_radius=78;
863  // double tpc_halfz=
864 
865  double ifc_radius = 83.5 * cm;
866  double ofc_radius = 254.5 * cm;
867  double tpc_halfz = 250 * cm;
868 
869  aliceModel = new AnalyticFieldModel(ifc_radius / cm, ofc_radius / cm, tpc_halfz / cm, scalefactor);
870  double totalcharge = 0;
871  double localcharge = 0;
872 
873  TVector3 pos;
874  for (int ifr = 0; ifr < nr; ifr++)
875  {
876  for (int ifphi = 0; ifphi < nphi; ifphi++)
877  {
878  for (int ifz = 0; ifz < nz; ifz++)
879  {
880  pos = GetCellCenter(ifr, ifphi, ifz);
881  pos = pos * (1.0 / cm); // divide by cm so we're in units of cm when we query the charge model.
882  double vol = step.Z() * step.Phi() * (2 * (ifr * step.Perp() + rmin) + step.Perp()) * step.Perp();
883  // if(debugFlag()) printf("%d: AnnularFieldSim::load_analytic_spacecharge adding Q=%f into cell (%d,%d,%d)\n",__LINE__,qbin,i,j,k,localr,localphi,localz);
884  localcharge = vol * aliceModel->Rho(pos); // TODO: figure out what units this is in.
885  totalcharge += localcharge;
886  // q->Add(ifr, ifphi, ifz, localcharge); //scalefactor must be applied to charge _and_ field, and so is handled in the aliceModel code.
887  q->AddChargeInBin(ifr, ifphi, ifz, localcharge); // scalefactor must be applied to charge _and_ field, and so is handled in the aliceModel code.
888  }
889  }
890  }
891  printf("AnnularFieldSim::load_analytic_spacecharge: Total charge Q=%E Coulombs\n", totalcharge);
892 
893  if (lookupCase == HybridRes)
894  {
895  // go through the q array and build q_lowres.
896  for (int i = 0; i < q_lowres->Length(); i++)
897  {
898  *(q_lowres->GetFlat(i)) = 0;
899  }
900 
901  // fill our low-res
902  // note that this assumes the last bin is short or normal length, not long.
903  for (int ifr = 0; ifr < nr; ifr++)
904  {
905  int r_low = ifr / r_spacing; // index of our l-bin is just the integer division of the index of our f-bin
906  for (int ifphi = 0; ifphi < nphi; ifphi++)
907  {
908  int phi_low = ifphi / phi_spacing;
909  for (int ifz = 0; ifz < nz; ifz++)
910  {
911  int z_low = ifz / z_spacing;
912  q_lowres->Add(r_low, phi_low, z_low, q->GetChargeInBin(ifr, ifphi, ifz));
913  }
914  }
915  }
916  }
917  return;
918 }
919 
920 void AnnularFieldSim::loadEfield(const std::string &filename, const std::string &treename, int zsign)
921 {
922  // prep variables so that loadField can just iterate over the tree entries and fill our selected tree agnostically
923  // assumes file stores fields as V/cm.
924  TFile fieldFile(filename.c_str(), "READ");
925  TTree *fTree;
926  fieldFile.GetObject(treename.c_str(), fTree);
927  Efieldname = "E:" + filename + ":" + treename;
928  // sprintf(Efieldname,"E:%s:%s",filename,treename);
929  float r, z; // coordinates
930  float fr, fz, fphi; // field components at that coordinate
931  fTree->SetBranchAddress("r", &r);
932  fTree->SetBranchAddress("er", &fr);
933  fTree->SetBranchAddress("z", &z);
934  fTree->SetBranchAddress("ez", &fz);
935  // phi would go here if we had it.
936  fphi = 0; // no phi components yet.
937  // float phi = 0.;
938  // phi += 1;
939  // phi = 0; //satisfy picky racf compiler
940  loadField(&Eexternal, fTree, &r, nullptr, &z, &fr, &fphi, &fz, V / cm, zsign);
941  fieldFile.Close();
942  return;
943 }
945 {
946  // prep variables so that loadField can just iterate over the tree entries and fill our selected tree agnostically
947  // assumes file stores field as Tesla.
948  TFile fieldFile(filename.c_str(), "READ");
949  TTree *fTree;
950  fieldFile.GetObject(treename.c_str(), fTree);
951  Bfieldname = "B:" + filename + ":" + treename;
952  float r, z; // coordinates
953  float fr, fz, fphi; // field components at that coordinate
954  fTree->SetBranchAddress("r", &r);
955  fTree->SetBranchAddress("br", &fr);
956  fTree->SetBranchAddress("z", &z);
957  fTree->SetBranchAddress("bz", &fz);
958  // phi would go here if we had it.
959  fphi = 0; // no phi components yet.
960  // float phi = 0.;
961  // phi += 1;
962  // phi = 0; //satisfy picky racf compiler
963  loadField(&Bfield, fTree, &r, nullptr, &z, &fr, &fphi, &fz, Tesla, 1);
964  fieldFile.Close();
965 
966  return;
967 }
968 
969 void AnnularFieldSim::load3dBfield(const std::string &filename, const std::string &treename, int zsign, float scale)
970 {
971  // prep variables so that loadField can just iterate over the tree entries and fill our selected tree agnostically
972  // assumes file stores field as Tesla.
973  TFile fieldFile(filename.c_str(), "READ");
974  TTree *fTree;
975  fieldFile.GetObject(treename.c_str(), fTree);
976  Bfieldname = "B(3D):" + filename + ":" + treename + " Scale =" + boost::str(boost::format("%2.2f") % scale);
977  // sprintf(Bfieldname,"B(3D):%s:%s Scale=%2.2f",filename,treename,scale);
978  float r, z, phi; // coordinates
979  float fr, fz, fphi; // field components at that coordinate
980  fTree->SetBranchAddress("rho", &r);
981  fTree->SetBranchAddress("brho", &fr);
982  fTree->SetBranchAddress("z", &z);
983  fTree->SetBranchAddress("bz", &fz);
984  fTree->SetBranchAddress("phi", &phi);
985  fTree->SetBranchAddress("bphi", &fphi);
986  loadField(&Bfield, fTree, &r, nullptr, &z, &fr, &fphi, &fz, Tesla * scale, zsign);
987  return;
988 }
989 
990 void AnnularFieldSim::loadField(MultiArray<TVector3> **field, TTree *source, float *rptr, float *phiptr, float *zptr, float *frptr, float *fphiptr, float *fzptr, float fieldunit, int zsign)
991 {
992  // we're loading a tree of unknown size and spacing -- and possibly uneven spacing -- into our local data.
993  // formally, we might want to interpolate or otherwise weight, but for now, carve this into our usual bins, and average, similar to the way we load spacecharge.
994 
995  bool phiSymmetry = (phiptr == nullptr); // if the phi pointer is zero, assume phi symmetry.
996  int lowres_factor = 10; // to fill in gaps, we group together loweres^3 cells into one block and use that average.
997  printf("loading field from %f<z<%f\n", zmin, zmax);
998  TH3F *htEntries = new TH3F("htentries", "num of entries in the field loading", nphi, 0, M_PI * 2.0, nr, rmin, rmax, nz, zmin, zmax);
999  TH3F *htSum[3];
1000  TH3F *htEntriesLow = new TH3F("htentrieslow", "num of lowres entries in the field loading", nphi / lowres_factor + 1, 0, M_PI * 2.0, nr / lowres_factor + 1, rmin, rmax, nz / lowres_factor + 1, zmin, zmax);
1001  TH3F *htSumLow[3];
1002  char axis[] = "rpz";
1003  for (int i = 0; i < 3; i++)
1004  {
1005  htSum[i] = new TH3F(Form("htsum%d", i), Form("sum of %c-axis entries in the field loading", *(axis + i)), nphi, 0, M_PI * 2.0, nr, rmin, rmax, nz, zmin, zmax);
1006  htSumLow[i] = new TH3F(Form("htsumlow%d", i), Form("sum of low %c-axis entries in the field loading", *(axis + i)), nphi / lowres_factor + 1, 0, M_PI * 2.0, nr / lowres_factor + 1, rmin, rmax, nz / lowres_factor + 1, zmin, zmax);
1007  }
1008 
1009  int nEntries = source->GetEntries();
1010  for (int i = 0; i < nEntries; i++)
1011  { // could probably do this with an iterator
1012  source->GetEntry(i);
1013  float zval = *zptr * zsign; // right now, need the ability to flip the sign of the z coordinate.
1014  // note that the z sign also needs to affect the field sign in that direction, which is handled outside in the z components of the fills
1015  // if we aren't asking for phi symmetry, build just the one phi strip
1016  if (!phiSymmetry)
1017  {
1018  assert(phiptr);
1019  htEntries->Fill(*phiptr, *rptr, zval); // for legacy reasons this histogram, like others, goes phi-r-z.
1020  htSum[0]->Fill(*phiptr, *rptr, zval, *frptr * fieldunit);
1021  htSum[1]->Fill(*phiptr, *rptr, zval, *fphiptr * fieldunit);
1022  htSum[2]->Fill(*phiptr, *rptr, zval, *fzptr * fieldunit * zsign);
1023  htEntriesLow->Fill(*phiptr, *rptr, zval); // for legacy reasons this histogram, like others, goes phi-r-z.
1024  htSumLow[0]->Fill(*phiptr, *rptr, zval, *frptr * fieldunit);
1025  htSumLow[1]->Fill(*phiptr, *rptr, zval, *fphiptr * fieldunit);
1026  htSumLow[2]->Fill(*phiptr, *rptr, zval, *fzptr * fieldunit * zsign);
1027  }
1028  else
1029  { // if we do have phi symmetry, build every phi strip using this one.
1030  for (int j = 0; j < nphi; j++)
1031  {
1032  htEntries->Fill(j * step.Phi(), *rptr, zval); // for legacy reasons this histogram, like others, goes phi-r-z.
1033  htSum[0]->Fill(j * step.Phi(), *rptr, zval, *frptr * fieldunit);
1034  htSum[1]->Fill(j * step.Phi(), *rptr, zval, *fphiptr * fieldunit);
1035  htSum[2]->Fill(j * step.Phi(), *rptr, zval, *fzptr * fieldunit * zsign);
1036  htEntriesLow->Fill(j * step.Phi(), *rptr, zval); // for legacy reasons this histogram, like others, goes phi-r-z.
1037  htSumLow[0]->Fill(j * step.Phi(), *rptr, zval, *frptr * fieldunit);
1038  htSumLow[1]->Fill(j * step.Phi(), *rptr, zval, *fphiptr * fieldunit);
1039  htSumLow[2]->Fill(j * step.Phi(), *rptr, zval, *fzptr * fieldunit * zsign);
1040  }
1041  }
1042  }
1043  // now we just divide and fill our local plots (which should eventually be stored as histograms, probably) with the values from each hist cell:
1044  int nemptybins = 0;
1045  for (int i = 0; i < nphi; i++)
1046  {
1047  for (int j = 0; j < nr; j++)
1048  {
1049  for (int k = 0; k < nz; k++)
1050  {
1051  TVector3 cellcenter = GetCellCenter(j, i, k);
1052  int bin = htEntries->FindBin(FilterPhiPos(cellcenter.Phi()), cellcenter.Perp(), cellcenter.Z());
1053  TVector3 fieldvec(htSum[0]->GetBinContent(bin), htSum[1]->GetBinContent(bin), htSum[2]->GetBinContent(bin));
1054  fieldvec = fieldvec * (1.0 / htEntries->GetBinContent(bin));
1055  if (htEntries->GetBinContent(bin) < 0.99)
1056  {
1057  // no entries here!
1058  nemptybins++;
1059  }
1060  // have to rotate this to the proper direction.
1061  fieldvec.RotateZ(FilterPhiPos(cellcenter.Phi())); // rcc caution. Does this rotation shift the sense of 'up'?
1062  (*field)->Set(j, i, k, fieldvec);
1063  }
1064  }
1065  }
1066  if (nemptybins > 0)
1067  {
1068  printf("found %d empty bins when constructing %s. Filling with lower resolution.\n", nemptybins, (*field == Bfield) ? "Bfield" : "Eexternal");
1069  for (int i = 0; i < nphi; i++)
1070  {
1071  for (int j = 0; j < nr; j++)
1072  {
1073  for (int k = 0; k < nz; k++)
1074  {
1075  TVector3 cellcenter = GetCellCenter(j, i, k);
1076  int bin = htEntries->FindBin(FilterPhiPos(cellcenter.Phi()), cellcenter.Perp(), cellcenter.Z());
1077  if (htEntries->GetBinContent(bin) == 0)
1078  {
1079  int lowbin = htEntriesLow->FindBin(FilterPhiPos(cellcenter.Phi()), cellcenter.Perp(), cellcenter.Z());
1080  TVector3 fieldvec(htSumLow[0]->GetBinContent(lowbin), htSumLow[1]->GetBinContent(lowbin), htSumLow[2]->GetBinContent(lowbin));
1081  fieldvec = fieldvec * (1.0 / htEntriesLow->GetBinContent(lowbin));
1082  if (htEntriesLow->GetBinContent(lowbin) < 0.99)
1083  {
1084  printf("not enough entries in source to fill fieldmap. None near r=%f, phi=%f, z=%f. Pick lower granularity!\n",
1085  cellcenter.Perp(), FilterPhiPos(cellcenter.Phi()), cellcenter.Z());
1086  assert(1 == 2);
1087  }
1088  // have to rotate this to the proper direction.
1089  fieldvec.RotateZ(FilterPhiPos(cellcenter.Phi())); // rcc caution. Does this rotation shift the sense of 'up'?
1090  (*field)->Set(j, i, k, fieldvec);
1091  }
1092  }
1093  }
1094  }
1095  }
1096  return;
1097 }
1098 
1099 void AnnularFieldSim::load_spacecharge(const std::string &filename, const std::string &histname, float zoffset, float chargescale, float cmscale, bool isChargeDensity)
1100 {
1101  TFile *f = TFile::Open(filename.c_str());
1102  TH3 *scmap;
1103  f->GetObject(histname.c_str(), scmap);
1104  std::cout << "Loading spacecharge from '" << filename
1105  << "'. Seeking histname '" << histname << "'" << std::endl;
1106  chargesourcename = filename + ":" + histname;
1107  // sprintf(chargesourcename,"%s:%s",filename,histname);
1108  load_spacecharge(scmap, zoffset, chargescale, cmscale, isChargeDensity, chargesourcename.c_str());
1109  f->Close();
1110  return;
1111 }
1112 
1113 void AnnularFieldSim::load_and_resample_spacecharge(int new_nphi, int new_nr, int new_nz, const std::string &filename, const std::string &histname, float zoffset, float chargescale, float cmscale, bool isChargeDensity)
1114 {
1115  TFile *f = TFile::Open(filename.c_str());
1116  TH3 *scmap;
1117  f->GetObject(histname.c_str(), scmap);
1118  std::cout << "Resampling spacecharge from '" << filename
1119  << "'. Seeking histname '" << histname << "'" << std::endl;
1120  chargesourcename = filename + ":" + histname;
1121  load_and_resample_spacecharge(new_nphi, new_nr, new_nz, scmap, zoffset, chargescale, cmscale, isChargeDensity);
1122  f->Close();
1123  return;
1124 }
1125 
1126 void AnnularFieldSim::load_and_resample_spacecharge(int new_nphi, int new_nr, int new_nz, TH3 *hist, float zoffset, float chargescale, float cmscale, bool isChargeDensity)
1127 {
1128  // load spacecharge densities from a histogram, where scalefactor translates into local units of C/cm^3
1129  // and cmscale translate (hist coord) --> (hist position in cm)
1130  // noting that the histogram limits may differ from the simulation size, and have different granularity
1131  // hist is assumed/required to be x=phi, y=r, z=z
1132  // z offset 'drifts' the charge by that distance toward z=0.
1133 
1134  // Get dimensions of input
1135  float hrmin = hist->GetYaxis()->GetXmin();
1136  float hrmax = hist->GetYaxis()->GetXmax();
1137  float hphimin = hist->GetXaxis()->GetXmin();
1138  float hphimax = hist->GetXaxis()->GetXmax();
1139  float hzmin = hist->GetZaxis()->GetXmin();
1140  float hzmax = hist->GetZaxis()->GetXmax();
1141 
1142  // Get number of bins in each dimension
1143  int hrn = hist->GetNbinsY();
1144  int hphin = hist->GetNbinsX();
1145  int hzn = hist->GetNbinsZ();
1146  printf("From histogram, native r dimensions: %f<r<%f, hrn (Y axis)=%d\n", hrmin, hrmax, hrn);
1147  printf("From histogram, native phi dimensions: %f<phi<%f, hrn (X axis)=%d\n", hphimin, hphimax, hphin);
1148  printf("From histogram, native z dimensions: %f<z<%f, hrn (Z axis)=%d\n", hzmin, hzmax, hzn);
1149 
1150  // do some computation of steps:
1151  float hrstep = (hrmax - hrmin) / hrn;
1152  float hphistep = (hphimax - hphimin) / hphin;
1153  float hzstep = (hzmax - hzmin) / hzn;
1154  float halfrstep = hrstep * 0.5;
1155  float halfphistep = hphistep * 0.5;
1156  float halfzstep = hzstep * 0.5;
1157 
1158  // All we have to do here is resample it so it fits into our expected grid, then pass it to the loader.
1159 
1160  // 1) convert the existing histogram to density if it isn't already:
1161  if (!isChargeDensity)
1162  {
1163  for (int i = 0; i < hphin; i++)
1164  {
1165  float phi = hphimin + hphistep * i;
1166  for (int j = 0; j < hrn; j++)
1167  {
1168  float r = hrmin + hrstep * j;
1169  for (int k = 0; k < hzn; k++)
1170  {
1171  float z = hzmin + hzstep * k;
1172  int bin = hist->FindBin(phi + halfphistep, r + halfrstep, z + halfzstep);
1173  double vol = hzstep * hphistep * (r + halfrstep) * hrstep;
1174  hist->SetBinContent(bin, hist->GetBinContent(bin) / vol);
1175  }
1176  }
1177  }
1178  }
1179  // I ought to consider a subtler approach that better preserves the overall integral, rather than just skipping the interpolation of the outer edges:
1180  TH3F *resampled = new TH3F("resampled", "resampled charge", new_nphi, hphimin, hphimax, new_nr, hrmin, hrmax, new_nz, hzmin, hzmax);
1181  float new_phistep = (hphimax - hphimin) / new_nphi;
1182  float new_rstep = (hrmax - hrmin) / new_nr;
1183  float new_zstep = (hzmax - hzmin) / new_nz;
1184 
1185  // 2 resample with interpolation on the interior, and with nearest-bin for the edge bins
1186  for (int i = 0; i < new_nphi; i++)
1187  {
1188  float phi = hphimin + new_phistep * (i + 0.5); // bin center of resampled
1189  float hphi = (phi - hphimin) / hphistep; // coord in source hist
1190  bool phisafe = ((hphi - 0.75) > 0) && ((hphi + 0.75) < hphin); // we're past the halfway point of the lowest bin, and short of the halfway point of the highest bin.
1191  for (int j = 0; j < new_nr; j++)
1192  {
1193  float r = hrmin + new_rstep * (j + 0.5); // bin center of resampled
1194  float hr = (r - hrmin) / hrstep; // coord in source hist
1195  bool rsafe = ((hr - 0.5) > 0) && ((hr + 0.5) < hrn); // we're past the halfway point of the lowest bin, and short of the halfway point of the highest bin.
1196  for (int k = 0; k < new_nz; k++)
1197  {
1198  float z = hzmin + new_zstep * (k + 0.5); // bin center of resampled
1199  float hz = (z - hzmin) / hzstep; // coord in source hist
1200  bool zsafe = ((hz - 0.5) > 0) && ((hz + 0.5) < hzn); // we're past the halfway point of the lowest bin, and short of the halfway point of the highest bin.
1201  // check if we need to do a bin lookup:
1202  if (phisafe && rsafe && zsafe)
1203  {
1204  // printf("resampling (%d,%d,%d) from (%2.2f/%d,%2.2f/%d,%2.2f/%d) p:(%1.1f<%1.1f<%1.1f) r:(%1.1f<%1.1f<%1.1f) z:(%1.1f<%1.1f<%1.1f)\n",i,j,k,hphi,hphin,hr,hrn,hz,hzn,hphimin,phi,hphimax,hrmin,r,hrmax,hzmin,z,hzmax);
1205  resampled->Fill(phi, r, z, hist->Interpolate(phi, r, z)); // leave as a density
1206  }
1207  else
1208  {
1209  int bin = hist->FindBin(phi, r, z);
1210  resampled->Fill(phi, r, z, hist->GetBinContent(bin));
1211  }
1212  } // k (z loop)
1213  } // j (r loop)
1214  } // i (phi loop)
1215  printf("here in load_and_resample_spacecharge, I am going to call load_spacecharge.\n");
1216  load_spacecharge(resampled, zoffset, chargescale, cmscale, true);
1217 }
1218 
1219 void AnnularFieldSim::load_spacecharge(TH3 *hist, float zoffset, float chargescale, float cmscale, bool isChargeDensity, const char *inputchargestring)
1220 {
1221  // new plan: use ChargeMapReader:
1222  if (abs(zoffset) > 0.001)
1223  {
1224  printf("nonzero zoffset given (%E) but new spacecharge loader can't deal with that. Failing.\n", zoffset);
1225  assert(false);
1226  }
1227  if (isChargeDensity)
1228  {
1229  printf("Input dataset is flagged as recording density not total charge, but new loader can't deal with that Failing..\n");
1230  assert(false);
1231  }
1232  q->ReadSourceCharge(hist, cmscale, chargescale);
1233 
1234  sprintf(chargestring, "SC loaded externally: %s.", inputchargestring);
1235  return;
1236 }
1237 
1238 void AnnularFieldSim::load_digital_current(TH3 *hist, TH2 *gainHist, float chargescale, float cmscale, const char *inputchargestring)
1239 {
1240  q->ReadSourceAdc(hist, gainHist, cmscale, chargescale);
1241 
1242  sprintf(chargestring, "SC loaded externally: %s.", inputchargestring);
1243  return;
1244 }
1245 
1247 {
1248  // save a histogram giving the spacecharge in our local coordinates, as we'll be using them..
1249  TH3F *hsc = new TH3F("hInternalSpaceCharge", "Internal Charge Histogram;phi(rad);r(cm);z(cm)", nphi, 0, phispan, nr, rmin, rmax, nz, zmin, zmax);
1250  for (int i = 0; i < nphi; i++)
1251  {
1252  float phi = 0 + step.Phi() * (i + 0.5);
1253  for (int j = 0; j < nr; j++)
1254  {
1255  float r = rmin + step.Perp() * (j + 0.5);
1256  for (int k = 0; k < nz; k++)
1257  {
1258  float z = zmin + step.Z() * (k + 0.5);
1259  hsc->Fill(phi, r, z, q->GetChargeAtPosition(r, phi, z));
1260  // old version: hsc->Fill(phi,r,z,q->Get(j,i,k));
1261  }
1262  }
1263  }
1264  hsc->SaveAs(filename.c_str());
1265  return;
1266 }
1267 
1268 void AnnularFieldSim::add_testcharge(float r, float phi, float z, float coulombs)
1269 {
1270  q->AddChargeAtPosition(r, phi, z, coulombs * C);
1271  return;
1272  /*
1273  int rcell, phicell, zcell;
1274 
1275  //translate to which cell we're in:
1276  BoundsCase rb = GetRindexAndCheckBounds(r, &rcell);
1277  BoundsCase pb = GetPhiIndexAndCheckBounds(phi, &phicell);
1278  BoundsCase zb = GetZindexAndCheckBounds(z, &zcell);
1279  //really, we would like to confirm that the cell we find is in the charge map region, not the field map region.
1280  if (rb == BoundsCase::OutOfBounds || pb == BoundsCase::OutOfBounds || zb == BoundsCase::OutOfBounds)
1281  {
1282  printf("placing %f Coulombs at (r,p,z)=(%f,%f,%f). Caution that that is out of the roi.\n", coulombs, r, phi, z);
1283  }
1284  if (rcell < 0 || phicell < 0 || zcell < 0)
1285  {
1286  printf("Tried placing %f Coulombs at (r,p,z)=(%f,%f,%f). That is outside of the charge map region.\n", coulombs, r, phi, z);
1287  return;
1288  }
1289 
1290  q->Add(rcell, phicell, zcell, coulombs * C);
1291  if (lookupCase == HybridRes)
1292  {
1293  printf("write the code you didn't write earlier about reloading the lowres map.\n");
1294  assert(1 == 2);
1295  }
1296 
1297  return;
1298  */
1299 }
1300 
1301 /*
1302 void AnnularFieldSim::populate_analytic_fieldmap(){
1303  //sum the E field at every point in the region of interest
1304  // remember that Efield uses relative indices
1305  printf("in pop_analytic_fieldmap, n=(%d,%d,%d)\n",nr,nphi,nz);
1306 
1307  TVector3 localF;//holder for the summed field at the current position.
1308  for (int ir=rmin_roi;ir<rmax_roi;ir++){
1309  for (int iphi=phimin_roi;iphi<phimax_roi;iphi++){
1310  for (int iz=zmin_roi;iz<zmax_roi;iz++){
1311  localF=aliceModel->E(GetCellCenter(ir, iphi, iz))+Eexternal->Get(ir-rmin_roi,iphi-phimin_roi,iz-zmin_roi);
1312  Efield->Set(ir-rmin_roi,iphi-phimin_roi,iz-zmin_roi,localF);
1313  //if (localF.Mag()>1e-9)
1314  if(debugFlag()) printf("%d: AnnularFieldSim::populate_analytic_fieldmap fieldmap@ (%d,%d,%d) mag=%f\n",__LINE__,ir,iphi,iz,localF.Mag());
1315  }
1316  }
1317  }
1318  return;
1319 }
1320 */
1321 
1323 {
1324  // sum the E field at every point in the region of interest
1325  // remember that Efield uses relative indices
1326  printf("in pop_fieldmap, n=(%d,%d,%d)\n", nr, nphi, nz);
1327 
1328  printf("populating fieldmap for (%dx%dx%d) grid with (%dx%dx%d) source \n", nr_roi, nphi_roi, nz_roi, nr, nphi, nz);
1329  if (truncation_length > 0)
1330  {
1331  printf(" ==> truncating anything more than %d cells away\n", truncation_length);
1332  }
1333  unsigned long long totalelements = nr_roi;
1334  totalelements *= nphi_roi;
1335  totalelements *= nz_roi; // breaking up this multiplication prevents a 32bit math overflow
1336  unsigned long long percent = totalelements / 100 * debug_npercent;
1337  printf("total elements = %llu\n", totalelements * nr * nphi * nz);
1338 
1339  int el = 0;
1340 
1341  TVector3 localF; // holder for the summed field at the current position.
1342  for (int ir = rmin_roi; ir < rmax_roi; ir++)
1343  {
1344  for (int iphi = phimin_roi; iphi < phimax_roi; iphi++)
1345  {
1346  for (int iz = zmin_roi; iz < zmax_roi; iz++)
1347  {
1348  localF = sum_field_at(ir, iphi, iz); // asks in global coordinates
1349  if (!(el % percent))
1350  {
1351  printf("populate_fieldmap %d%%: ", (int) (debug_npercent * el / percent));
1352  printf("sum_field_at (ir=%d,iphi=%d,iz=%d) gives (%E,%E,%E)\n",
1353  ir, iphi, iz, localF.X(), localF.Y(), localF.Z());
1354  }
1355  el++;
1356 
1357  Efield->Set(ir - rmin_roi, iphi - phimin_roi, iz - zmin_roi, localF); // sets in roi coordinates.
1358  // if (localF.Mag()>1e-9)
1359  // if(debugFlag()) printf("%d: AnnularFieldSim::populate_fieldmap fieldmap@ (%d,%d,%d) mag=%f\n",__LINE__,ir,iphi,iz,localF.Mag());
1360  }
1361  }
1362  }
1363  return;
1364 }
1365 
1367 {
1368  // with 'f' being the position the field is being measured at, and 'o' being the position of the charge generating the field.
1369  // remember the 'f' part of Epartial uses relative indices.
1370  // TVector3 (*f)[fx][fy][fz][ox][oy][oz]=field_;
1371  // printf("populating lookup for (%dx%dx%d)x(%dx%dx%d) grid\n",fx,fy,fz,ox,oy,oz);
1372 
1373  if (lookupCase == Full3D)
1374  {
1375  printf("lookupCase==Full3D\n");
1376 
1378  }
1379  else if (lookupCase == HybridRes)
1380  {
1381  printf("lookupCase==HybridRes\n");
1384  }
1385  else if (lookupCase == PhiSlice)
1386  {
1387  printf("Populating lookup: lookupCase==PhiSlice\n");
1389  }
1390  else if (lookupCase == Analytic)
1391  {
1392  printf("Populating lookup: lookupCase==Analytic ===> skipping!\n");
1393  }
1394  else if (lookupCase == NoLookup)
1395  {
1396  printf("Populating lookup: lookupCase==NoLookup ===> skipping!\n");
1397  }
1398  else
1399  {
1400  assert(1 == 2);
1401  }
1402  return;
1403 }
1404 
1406 {
1407  // with 'f' being the position the field is being measured at, and 'o' being the position of the charge generating the field.
1408  // remember the 'f' part of Epartial uses relative indices.
1409  // TVector3 (*f)[fx][fy][fz][ox][oy][oz]=field_;
1410  printf("populating full lookup table for (%dx%dx%d)x(%dx%dx%d) grid\n",
1412  unsigned long long totalelements = (rmax_roi - rmin_roi);
1413  totalelements *= (phimax_roi - phimin_roi);
1414  totalelements *= (zmax_roi - zmin_roi);
1415  totalelements *= nr;
1416  totalelements *= nphi;
1417  totalelements *= nz; // breaking up this multiplication prevents a 32bit math overflow
1418  unsigned long long percent = totalelements / 100 * debug_npercent;
1419  printf("total elements = %llu\n", totalelements);
1420  TVector3 at(1, 0, 0);
1421  TVector3 from(1, 0, 0);
1422  TVector3 zero(0, 0, 0);
1423 
1424  int el = 0;
1425  for (int ifr = rmin_roi; ifr < rmax_roi; ifr++)
1426  {
1427  for (int ifphi = phimin_roi; ifphi < phimax_roi; ifphi++)
1428  {
1429  for (int ifz = zmin_roi; ifz < zmax_roi; ifz++)
1430  {
1431  at = GetCellCenter(ifr, ifphi, ifz);
1432  for (int ior = 0; ior < nr; ior++)
1433  {
1434  for (int iophi = 0; iophi < nphi; iophi++)
1435  {
1436  for (int ioz = 0; ioz < nz; ioz++)
1437  {
1438  el++;
1439  if (!(el % percent))
1440  {
1441  printf("populate_full3d_lookup %d%%\n", (int) (debug_npercent * el / percent));
1442  }
1443  from = GetCellCenter(ior, iophi, ioz);
1444 
1445  //*f[ifx][ify][ifz][iox][ioy][ioz]=cacl_unit_field(at,from);
1446  // printf("calc_unit_field...\n");
1447  if (ifr == ior && ifphi == iophi && ifz == ioz)
1448  {
1449  Epartial->Set(ifr - rmin_roi, ifphi - phimin_roi, ifz - zmin_roi, ior, iophi, ioz, zero);
1450  }
1451  else
1452  {
1453  Epartial->Set(ifr - rmin_roi, ifphi - phimin_roi, ifz - zmin_roi, ior, iophi, ioz, calc_unit_field(at, from));
1454  }
1455  }
1456  }
1457  }
1458  }
1459  }
1460  }
1461  return;
1462 }
1463 
1465 {
1466  TVector3 at(1, 0, 0);
1467  TVector3 from(1, 0, 0);
1468  TVector3 zero(0, 0, 0);
1469 
1470  // populate_highres_lookup();
1471  int r_highres_dist = (nr_high - 1) / 2;
1472  int phi_highres_dist = (nphi_high - 1) / 2;
1473  int z_highres_dist = (nz_high - 1) / 2;
1474 
1475  // number of fbins contained in the 26 weirdly-shaped edge regions (and one center region which we won't use)
1476  static int nfbinsin[3][3][3];
1477  for (auto &i : nfbinsin)
1478  {
1479  for (int j = 0; j < 3; j++)
1480  {
1481  for (int k = 0; k < 3; k++)
1482  {
1483  i[j][k] = 0; // we could count total volume, but without knowing the charge prior, it's not clear that'd be /better/
1484  }
1485  }
1486  }
1487 
1488  // todo: if this runs too slowly, I can do geometry instead of looping over all the cells that are possibly in range
1489 
1490  // loop over all the f-bins in the roi:
1491  TVector3 currentf, newf; // the averaged field vector without, and then with the new contribution from the f-bin being considered.
1492  for (int ifr = rmin_roi; ifr < rmax_roi; ifr++)
1493  {
1494  int r_parentlow = floor((ifr - r_highres_dist) / (r_spacing * 1.0)); // l-bin partly enclosed in our high-res region
1495  int r_parenthigh = floor((ifr + r_highres_dist) / (r_spacing * 1.0)) + 1; // definitely not enclosed in our high-res region
1496  int r_startpoint = r_parentlow * r_spacing; // the first f-bin of the lowest-r f-bin that our h-region touches. COuld be less than zero.
1497  int r_endpoint = r_parenthigh * r_spacing; // the first f-bin of the lowest-r l-bin after that that our h-region does not touch. could be larger than max.
1498 
1499  for (int ifphi = phimin_roi; ifphi < phimax_roi; ifphi++)
1500  {
1501  int phi_parentlow = floor(FilterPhiIndex(ifphi - phi_highres_dist) / (phi_spacing * 1.0)); // note this may have wrapped around
1502  bool phi_parentlow_wrapped = (ifphi - phi_highres_dist < 0);
1503  int phi_startpoint = phi_parentlow * phi_spacing; // the first f-bin of the lowest-z f-bin that our h-region touches.
1504  if (phi_parentlow_wrapped)
1505  {
1506  phi_startpoint -= nphi; // if we wrapped, re-wrap us so we're negative again
1507  }
1508 
1509  int phi_parenthigh = floor(FilterPhiIndex(ifphi + phi_highres_dist) / (phi_spacing * 1.0)) + 1; // note that this may have wrapped around
1510  bool phi_parenthigh_wrapped = (ifphi + phi_highres_dist >= nphi);
1511  int phi_endpoint = phi_parenthigh * phi_spacing;
1512  if (phi_parenthigh_wrapped)
1513  {
1514  phi_endpoint += nphi; // if we wrapped, re-wrap us so we're larger than nphi again. We use these relative coords to determine the position relative to the center of our h-region.
1515  }
1516 
1517  for (int ifz = zmin_roi; ifz < zmax_roi; ifz++)
1518  {
1519  // if(debugFlag()) printf("%d: AnnularFieldSim::populate_highres_lookup icell=(%d,%d,%d)\n",__LINE__,ifr,ifphi,ifz);
1520 
1521  int z_parentlow = floor((ifz - z_highres_dist) / (z_spacing * 1.0));
1522  int z_parenthigh = floor((ifz + z_highres_dist) / (z_spacing * 1.0)) + 1;
1523  int z_startpoint = z_parentlow * z_spacing; // the first f-bin of the lowest-z f-bin that our h-region touches.
1524  int z_endpoint = z_parenthigh * z_spacing; // the first f-bin of the lowest-z l-bin after that that our h-region does not touch.
1525 
1526  // our 'at' position, in global coords:
1527  at = GetCellCenter(ifr, ifphi, ifz);
1528  // define the farthest-away parent l-bin cells we're dealing with here:
1529  // note we're still in absolute coordinates
1530 
1531  // for every f-bin in the l-bins we're dealing with, figure out which relative highres bin it's in, and average the field vector into that bin's vector
1532  // note most of these relative bins have exactly one f-bin in them. It's only the edges that can get more.
1533  // note this is a running average: Anew=(Aold*Nold+V)/(Nold+1) and so on.
1534  // note also that we automatically skip f-bins that would've been out of the valid overall volume.
1535  for (int ir = r_startpoint; ir < r_endpoint; ir++)
1536  {
1537  // skip parts that are out of range:
1538  // could speed this up by moving this into the definition of start and endpoint.
1539  if (ir < 0)
1540  {
1541  ir = 0;
1542  }
1543  if (ir >= nr)
1544  {
1545  break;
1546  }
1547 
1548  int rbin = (ir - ifr) + r_highres_dist; // zeroth bin when we're at max distance below, etc.
1549  int rcell = 1;
1550  if (rbin <= 0)
1551  {
1552  rbin = 0;
1553  rcell = 0;
1554  }
1555  if (rbin >= nr_high)
1556  {
1557  rbin = nr_high - 1;
1558  rcell = 2;
1559  }
1560 
1561  for (int iphi = phi_startpoint; iphi < phi_endpoint; iphi++)
1562  {
1563  // no phi out-of-range checks since it's circular, but we provide ourselves a filtered version:
1564  int phiFilt = FilterPhiIndex(iphi);
1565  int phibin = (iphi - ifphi) + phi_highres_dist;
1566  int phicell = 1;
1567  if (phibin <= 0)
1568  {
1569  phibin = 0;
1570  phicell = 0;
1571  }
1572  if (phibin >= nphi_high)
1573  {
1574  phibin = nphi_high - 1;
1575  phicell = 2;
1576  }
1577  for (int iz = z_startpoint; iz < z_endpoint; iz++)
1578  {
1579  if (iz < 0)
1580  {
1581  iz = 0;
1582  }
1583  if (iz >= nz)
1584  {
1585  break;
1586  }
1587  int zbin = (iz - ifz) + z_highres_dist;
1588  int zcell = 1;
1589  if (zbin <= 0)
1590  {
1591  zbin = 0;
1592  zcell = 0;
1593  }
1594  if (zbin >= nz_high)
1595  {
1596  zbin = nz_high - 1;
1597  zcell = 2;
1598  }
1599  //'from' is in absolute coordinates
1600  from = GetCellCenter(ir, phiFilt, iz);
1601 
1602  nfbinsin[rcell][phicell][zcell]++;
1603  int nf = nfbinsin[rcell][phicell][zcell];
1604  // coordinates relative to the region of interest:
1605  int ir_rel = ifr - rmin_roi;
1606  int iphi_rel = ifphi - phimin_roi;
1607  int iz_rel = ifz - zmin_roi;
1608  if (zcell != 1 || rcell != 1 || phicell != 1)
1609  {
1610  // we're not in the center, so deal with our weird shapes by averaging:
1611  // but Epartial is in coordinates relative to the roi
1612  if (iphi_rel < 0)
1613  {
1614  printf("%d: Getting with phi=%d\n", __LINE__, iphi_rel);
1615  }
1616  currentf = Epartial_highres->Get(ir_rel, iphi_rel, iz_rel, rbin, phibin, zbin);
1617  // to keep this as the average, we multiply what's there back to its initial summed-but-not-divided value
1618  // then add our new value, and the divide the new sum by the total number of cells
1619  newf = (currentf * (nf - 1) + calc_unit_field(at, from)) * (1 / (nf * 1.0));
1620  Epartial_highres->Set(ir_rel, iphi_rel, iz_rel, rbin, phibin, zbin, newf);
1621  }
1622  else
1623  {
1624  // we're in the center cell, which means any f-bin that's not on the outer edge of our region:
1625  // calc_unit_field will return zero when at=from, so the center will be automatically zero here.
1626  if (ifr == rbin && ifphi == phibin && ifz == zbin)
1627  {
1628  Epartial_highres->Set(ir_rel, iphi_rel, iz_rel, rbin, phibin, zbin, zero);
1629  }
1630  else
1631  { // for extra carefulness, only calc the field if it's not self-to-self.
1632  newf = calc_unit_field(at, from);
1633  Epartial_highres->Set(ir_rel, iphi_rel, iz_rel, rbin, phibin, zbin, newf);
1634  }
1635  }
1636  }
1637  }
1638  }
1639  }
1640  }
1641  }
1642  return;
1643 }
1644 
1646 {
1647  TVector3 at(1, 0, 0);
1648  TVector3 from(1, 0, 0);
1649  TVector3 zero(0, 0, 0);
1650  int fphi_low, fphi_high, fz_low, fz_high; // edges of the outer l-bin
1651  int r_low, r_high, phi_low, phi_high, z_low, z_high; // edges of the inner l-bin
1652 
1653  // todo: add in handling if roi_low is wrap-around in phi
1654  for (int ifr = rmin_roi_low; ifr < rmax_roi_low; ifr++)
1655  {
1656  int fr_low = ifr * r_spacing;
1657  int fr_high = fr_low + r_spacing - 1;
1658  if (fr_high >= nr)
1659  {
1660  fr_high = nr - 1;
1661  }
1662  for (int ifphi = phimin_roi_low; ifphi < phimax_roi_low; ifphi++)
1663  {
1664  fphi_low = ifphi * phi_spacing;
1665  fphi_high = fphi_low + phi_spacing - 1;
1666  if (fphi_high >= nphi)
1667  {
1668  fphi_high = nphi - 1; // if our phi l-bins aren't evenly spaced, we need to catch that here.
1669  }
1670  for (int ifz = zmin_roi_low; ifz < zmax_roi_low; ifz++)
1671  {
1672  fz_low = ifz * z_spacing;
1673  fz_high = fz_low + z_spacing - 1;
1674  if (fz_high >= nz)
1675  {
1676  fz_high = nz - 1;
1677  }
1678  at = GetGroupCellCenter(fr_low, fr_high, fphi_low, fphi_high, fz_low, fz_high);
1679  // printf("ifr=%d, rlow=%d,rhigh=%d,r_spacing=%d\n",ifr,r_low,r_high,r_spacing);
1680  // if(debugFlag()) printf("%d: AnnularFieldSim::populate_lowres_lookup icell=(%d,%d,%d)\n",__LINE__,ifr,ifphi,ifz);
1681 
1682  for (int ior = 0; ior < nr_low; ior++)
1683  {
1684  r_low = ior * r_spacing;
1685  r_high = r_low + r_spacing - 1;
1686  int ir_rel = ifr - rmin_roi_low;
1687 
1688  if (r_high >= nr)
1689  {
1690  r_high = nr - 1;
1691  }
1692  for (int iophi = 0; iophi < nphi_low; iophi++)
1693  {
1694  phi_low = iophi * phi_spacing;
1695  phi_high = phi_low + phi_spacing - 1;
1696  if (phi_high >= nphi)
1697  {
1698  phi_high = nphi - 1;
1699  }
1700  int iphi_rel = ifphi - phimin_roi_low;
1701  for (int ioz = 0; ioz < nz_low; ioz++)
1702  {
1703  z_low = ioz * z_spacing;
1704  z_high = z_low + z_spacing - 1;
1705  if (z_high >= nz)
1706  {
1707  z_high = nz - 1;
1708  }
1709  int iz_rel = ifz - zmin_roi_low;
1710  from = GetGroupCellCenter(r_low, r_high, phi_low, phi_high, z_low, z_high);
1711 
1712  if (ifr == ior && ifphi == iophi && ifz == ioz)
1713  {
1714  Epartial_lowres->Set(ir_rel, iphi_rel, iz_rel, ior, iophi, ioz, zero);
1715  }
1716  else
1717  { // for extra carefulness, only calc the field if it's not self-to-self.
1718  Epartial_lowres->Set(ir_rel, iphi_rel, iz_rel, ior, iophi, ioz, calc_unit_field(at, from));
1719  }
1720 
1721  //*f[ifx][ify][ifz][iox][ioy][ioz]=cacl_unit_field(at,from);
1722  // printf("calc_unit_field...\n");
1723  // this calc's okay.
1724  }
1725  }
1726  }
1727  }
1728  }
1729  }
1730  return;
1731 }
1732 
1734 {
1735  // with 'f' being the position the field is being measured at, and 'o' being the position of the charge generating the field.
1736  // remember the 'f' part of Epartial uses relative indices.
1737  // TVector3 (*f)[fx][fy][fz][ox][oy][oz]=field_;
1738  printf("populating phislice lookup for (%dx%dx%d)x(%dx%dx%d) grid\n", nr_roi, 1, nz_roi, nr, nphi, nz);
1739  unsigned long long totalelements = nr; // nr*nphi*nz*nr_roi*nz_roi
1740  totalelements *= nphi;
1741  totalelements *= nz;
1742  totalelements *= nr_roi;
1743  totalelements *= nz_roi; // breaking up this multiplication prevents a 32bit math overflow
1744  unsigned long long percent = totalelements / 100 * debug_npercent;
1745  printf("total elements = %llu\n", totalelements);
1746  TVector3 at(1, 0, 0);
1747  TVector3 from(1, 0, 0);
1748  TVector3 zero(0, 0, 0);
1749 
1750  int el = 0;
1751  for (int ifr = rmin_roi; ifr < rmax_roi; ifr++)
1752  {
1753  for (int ifz = zmin_roi; ifz < zmax_roi; ifz++)
1754  {
1755  at = GetCellCenter(ifr, 0, ifz);
1756  for (int ior = 0; ior < nr; ior++)
1757  {
1758  for (int iophi = 0; iophi < nphi; iophi++)
1759  {
1760  for (int ioz = 0; ioz < nz; ioz++)
1761  {
1762  el++;
1763  from = GetCellCenter(ior, iophi, ioz);
1764  //*f[ifx][ify][ifz][iox][ioy][ioz]=cacl_unit_field(at,from);
1765  // printf("calc_unit_field...\n");
1766  if (ifr == ior && 0 == iophi && ifz == ioz)
1767  {
1768  if (!(el % percent))
1769  {
1770  printf("populate_phislice_lookup %d%%: ", (int) (debug_npercent * el / percent));
1771  printf("self-to-self is zero (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) gives (%E,%E,%E)\n",
1772  ior, iophi, ioz, ifr, ifz, zero.X(), zero.Y(), zero.Z());
1773  }
1774  Epartial_phislice->Set(ifr - rmin_roi, 0, ifz - zmin_roi, ior, iophi, ioz, zero);
1775  }
1776  else
1777  {
1778  TVector3 unitf = calc_unit_field(at, from);
1779  if (true)
1780  {
1781  if (!(el % percent))
1782  {
1783  printf("populate_phislice_lookup %d%%: ", (int) (debug_npercent * el / percent));
1784  printf("calc_unit_field (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) gives (%E,%E,%E)\n",
1785  ior, iophi, ioz, ifr, ifz, unitf.X(), unitf.Y(), unitf.Z());
1786  }
1787  }
1788 
1789  Epartial_phislice->Set(ifr - rmin_roi, 0, ifz - zmin_roi, ior, iophi, ioz, unitf); // the origin phi is relative to zero anyway.
1790  }
1791  }
1792  }
1793  }
1794  }
1795  }
1796  return;
1797 }
1798 
1799 void AnnularFieldSim::load_phislice_lookup(const char *sourcefile)
1800 {
1801  printf("loading phislice lookup for (%dx%dx%d)x(%dx%dx%d) grid from %s\n", nr_roi, 1, nz_roi, nr, nphi, nz, sourcefile);
1802  unsigned long long totalelements = nr; // nr*nphi*nz*nr_roi*nz_roi
1803  totalelements *= nphi;
1804  totalelements *= nz;
1805  totalelements *= nr_roi;
1806  totalelements *= nz_roi; // breaking up this multiplication prevents a 32bit math overflow
1807  unsigned long long percent = totalelements / 100 * debug_npercent;
1808  printf("total elements = %llu\n", totalelements);
1809 
1810  TFile *input = TFile::Open(sourcefile, "READ");
1811  TTree *tInfo;
1812  input->GetObject("info", tInfo);
1813  assert(tInfo);
1814 
1815  float file_rmin, file_rmax, file_zmin, file_zmax;
1816  int file_rmin_roi, file_rmax_roi, file_zmin_roi, file_zmax_roi;
1817  int file_nr, file_np, file_nz;
1818  tInfo->SetBranchAddress("rmin", &file_rmin);
1819  tInfo->SetBranchAddress("rmax", &file_rmax);
1820  tInfo->SetBranchAddress("zmin", &file_zmin);
1821  tInfo->SetBranchAddress("zmax", &file_zmax);
1822  tInfo->SetBranchAddress("rmin_roi_index", &file_rmin_roi);
1823  tInfo->SetBranchAddress("rmax_roi_index", &file_rmax_roi);
1824  tInfo->SetBranchAddress("zmin_roi_index", &file_zmin_roi);
1825  tInfo->SetBranchAddress("zmax_roi_index", &file_zmax_roi);
1826  tInfo->SetBranchAddress("nr", &file_nr);
1827  tInfo->SetBranchAddress("nphi", &file_np);
1828  tInfo->SetBranchAddress("nz", &file_nz);
1829  tInfo->GetEntry(0);
1830 
1831  printf("param\tobj\tfile\n");
1832  printf("nr\t%d\t%d\n", nr, file_nr);
1833  printf("np\t%d\t%d\n", nphi, file_np);
1834  printf("nz\t%d\t%d\n", nz, file_nz);
1835  printf("rmin\t%2.2f\t%2.2f\n", rmin, file_rmin);
1836  printf("rmax\t%2.2f\t%2.2f\n", rmax, file_rmax);
1837  printf("zmin\t%2.2f\t%2.2f\n", zmin, file_zmin);
1838  printf("zmax\t%2.2f\t%2.2f\n", zmax, file_zmax);
1839  printf("rmin_roi\t%d\t%d\n", rmin_roi, file_rmin_roi);
1840  printf("rmax_roi\t%d\t%d\n", rmax_roi, file_rmax_roi);
1841  printf("zmin_roi\t%d\t%d\n", zmin_roi, file_zmin_roi);
1842  printf("zmax_roi\t%d\t%d\n", zmax_roi, file_zmax_roi);
1843 
1844  if (file_rmin != rmin || file_rmax != rmax ||
1845  file_zmin != zmin || file_zmax != zmax ||
1846  file_rmin_roi != rmin_roi || file_rmax_roi != rmax_roi ||
1847  file_zmin_roi != zmin_roi || file_zmax_roi != zmax_roi ||
1848  file_nr != nr || file_np != nphi || file_nz != nz)
1849  {
1850  printf("file parameters do not match fieldsim parameters:\n");
1851 
1852  assert(1 == 4);
1853  }
1854 
1855  TTree *tLookup;
1856  input->GetObject("phislice", tLookup);
1857  assert(tLookup);
1858  int ior, ifr, iophi, ioz, ifz;
1859  TVector3 *unitf = nullptr;
1860  tLookup->SetBranchAddress("ir_source", &ior);
1861  tLookup->SetBranchAddress("ir_target", &ifr);
1862  tLookup->SetBranchAddress("ip_source", &iophi);
1863  // always zero: tLookup->SetBranchAddress("ip_target",&ifphi);
1864  tLookup->SetBranchAddress("iz_source", &ioz);
1865  tLookup->SetBranchAddress("iz_target", &ifz);
1866  tLookup->SetBranchAddress("Evec", &unitf); // assume field has units V/(C*cm)
1867 
1868  int el = 0;
1869  printf("%s has %lld entries\n", sourcefile, tLookup->GetEntries());
1870  for (int i = 0; i < (int) totalelements; i++)
1871  {
1872  el++;
1873  tLookup->GetEntry(i);
1874  // printf("loading i=%d\n",i);
1875  Epartial_phislice->Set(ifr - rmin_roi, 0, ifz - zmin_roi, ior, iophi, ioz, (*unitf) * (-1.0) * (V / (cm * C))); // load assuming field has units V/(C*cm), which is how we save it.
1876  // note that we save the gradient terms, not the field, hence we need to multiply by (-1.0)
1877  if (!(el % percent))
1878  {
1879  printf("load_phislice_lookup %d%%: ", (int) (debug_npercent * el / percent));
1880  printf("field from (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) is (%E,%E,%E)\n",
1881  ior, iophi, ioz, ifr, ifz, unitf->X(), unitf->Y(), unitf->Z());
1882  }
1883  }
1884 
1885  input->Close();
1886  return;
1887 }
1888 
1889 void AnnularFieldSim::save_phislice_lookup(const char *destfile)
1890 {
1891  printf("saving phislice lookup for (%dx%dx%d)x(%dx%dx%d) grid to %s\n", nr_roi, 1, nz_roi, nr, nphi, nz, destfile);
1892  unsigned long long totalelements = nr; // nr*nphi*nz*nr_roi*nz_roi
1893  totalelements *= nphi;
1894  totalelements *= nz;
1895  totalelements *= nr_roi;
1896  totalelements *= nz_roi; // breaking up this multiplication prevents a 32bit math overflow
1897  unsigned long long percent = totalelements / 100 * debug_npercent;
1898  printf("total elements = %llu\n", totalelements);
1899 
1900  TFile *output = TFile::Open(destfile, "RECREATE");
1901  output->cd();
1902 
1903  TTree *tInfo = new TTree("info", "Information about Lookup Table");
1904  tInfo->Branch("rmin", &rmin);
1905  tInfo->Branch("rmax", &rmax);
1906  tInfo->Branch("zmin", &zmin);
1907  tInfo->Branch("zmax", &zmax);
1908  tInfo->Branch("rmin_roi_index", &rmin_roi);
1909  tInfo->Branch("rmax_roi_index", &rmax_roi);
1910  tInfo->Branch("zmin_roi_index", &zmin_roi);
1911  tInfo->Branch("zmax_roi_index", &zmax_roi);
1912  tInfo->Branch("nr", &nr);
1913  tInfo->Branch("nphi", &nphi);
1914  tInfo->Branch("nz", &nz);
1915  tInfo->Fill();
1916  printf("info tree built.\n");
1917 
1918  TTree *tLookup = new TTree("phislice", "Phislice Lookup Table");
1919  int ior, ifr, iophi, ioz, ifz;
1920  TVector3 unitf;
1921  tLookup->Branch("ir_source", &ior);
1922  tLookup->Branch("ir_target", &ifr);
1923  tLookup->Branch("ip_source", &iophi);
1924  // always zero: tLookup->Branch("ip_target",&ifphi);
1925  tLookup->Branch("iz_source", &ioz);
1926  tLookup->Branch("iz_target", &ifz);
1927  tLookup->Branch("Evec", &unitf);
1928  printf("lookup tree built.\n");
1929 
1930  int el = 0;
1931  for (ifr = rmin_roi; ifr < rmax_roi; ifr++)
1932  {
1933  for (ifz = zmin_roi; ifz < zmax_roi; ifz++)
1934  {
1935  for (ior = 0; ior < nr; ior++)
1936  {
1937  for (iophi = 0; iophi < nphi; iophi++)
1938  {
1939  for (ioz = 0; ioz < nz; ioz++)
1940  {
1941  el++;
1942  unitf = Epartial_phislice->Get(ifr - rmin_roi, 0, ifz - zmin_roi, ior, iophi, ioz) * (-1 / (V / (C * cm))); // save in units of V/(C*cm) note that we introduce a -1 here for legcy reasons.
1943  if (true)
1944  {
1945  if (!(el % percent))
1946  {
1947  printf("save_phislice_lookup %d%%: ", (int) (debug_npercent * el / percent));
1948  printf("field from (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) is (%E,%E,%E)\n",
1949  ior, iophi, ioz, ifr, ifz, unitf.X(), unitf.Y(), unitf.Z());
1950  }
1951  }
1952 
1953  tLookup->Fill();
1954  }
1955  }
1956  }
1957  }
1958  }
1959  output->cd();
1960  tInfo->Write();
1961  tLookup->Write();
1962  // output->Write();
1963  output->Close();
1964  return;
1965 }
1966 
1967 void AnnularFieldSim::setFlatFields(float B, float E)
1968 {
1969  // these only cover the roi, but since we address them flat, we don't need to know that here.
1970  printf("AnnularFieldSim::setFlatFields(B=%f Tesla,E=%f V/cm)\n", B, E);
1971  printf("lengths: Eext=%d, Bfie=%d\n", Eexternal->Length(), Bfield->Length());
1972  char fieldstr[100];
1973  sprintf(fieldstr, "%f", E);
1974  Efieldname = "E:Flat:" + std::string(fieldstr);
1975  sprintf(fieldstr, "%f", B);
1976  Bfieldname = "B:Flat:" + std::string(fieldstr);
1977 
1978  Enominal = E * (V / cm);
1979  Bnominal = B * Tesla;
1980  for (int i = 0; i < Eexternal->Length(); i++)
1981  {
1982  Eexternal->GetFlat(i)->SetXYZ(0, 0, Enominal);
1983  }
1984  for (int i = 0; i < Bfield->Length(); i++)
1985  {
1986  Bfield->GetFlat(i)->SetXYZ(0, 0, Bnominal);
1987  }
1988  UpdateOmegaTau();
1989  return;
1990 }
1991 
1992 TVector3 AnnularFieldSim::sum_field_at(int r, int phi, int z)
1993 {
1994  // sum the E field over all nr by ny by nz cells of sources, at the global coordinate position r,phi,z.
1995  // note the specific position in Epartial is in relative coordinates.
1996  // if(debugFlag()) printf("%d: AnnularFieldSim::sum_field_at(r=%d,phi=%d, z=%d)\n",__LINE__,r,phi,z);
1997 
1998  /*
1999  for the near future:
2000  TVector3 sum=(sum_local_field_at(r, phi, z)
2001  + sum_nonlocal_field_at(r,phi,z)
2002  + Escale*Eexternal->Get(r-rmin_roi,phi-phimin_roi,z-zmin_roi));
2003  return sum;
2004  */
2005 
2006  TVector3 sum(0, 0, 0);
2007  if (lookupCase == Full3D)
2008  {
2009  sum += sum_full3d_field_at(r, phi, z);
2010  }
2011  else if (lookupCase == HybridRes)
2012  {
2013  sum += sum_local_field_at(r, phi, z);
2014  sum += sum_nonlocal_field_at(r, phi, z);
2015  }
2016  else if (lookupCase == PhiSlice)
2017  {
2018  sum += sum_phislice_field_at(r, phi, z);
2019  }
2020  else if (lookupCase == Analytic)
2021  {
2022  sum += aliceModel->E(GetCellCenter(r, phi, z));
2023  }
2024  else if (lookupCase == NoLookup)
2025  {
2026  // do nothing. We are forcibly assuming E from spacecharge is zero everywhere.
2027  }
2028  sum += Eexternal->Get(r - rmin_roi, phi - phimin_roi, z - zmin_roi);
2029  if (debugFlag())
2030  {
2031  printf("summed field at (%d,%d,%d)=(%f,%f,%f)\n", r, phi, z, sum.X(), sum.Y(), sum.Z());
2032  }
2033 
2034  return sum;
2035 }
2036 
2038 {
2039  // sum the E field over all nr by ny by nz cells of sources, at the specific position r,phi,z.
2040  // note the specific position in Epartial is in relative coordinates.
2041  // printf("AnnularFieldSim::sum_field_at(r=%d,phi=%d, z=%d)\n",r,phi,z);
2042  TVector3 sum(0, 0, 0);
2043  float rdist, phidist, zdist, remdist;
2044  for (int ir = 0; ir < nr; ir++)
2045  {
2046  if (truncation_length > 0)
2047  {
2048  rdist = abs(ir - r);
2049  remdist = sqrt(truncation_length * truncation_length - rdist * rdist);
2050  if (remdist < 0)
2051  {
2052  continue; // skip if we're too far away
2053  }
2054  }
2055  for (int iphi = 0; iphi < nphi; iphi++)
2056  {
2057  if (truncation_length > 0)
2058  {
2059  phidist = fmin(abs(iphi - phi), abs(abs(iphi - phi) - nphi)); // think about this in phi... rcc food for thought.
2060  remdist = sqrt(truncation_length * truncation_length - phidist * phidist);
2061  if (remdist < 0)
2062  {
2063  continue; // skip if we're too far away
2064  }
2065  }
2066  for (int iz = 0; iz < nz; iz++)
2067  {
2068  if (truncation_length > 0)
2069  {
2070  zdist = abs(iz - z);
2071  remdist = sqrt(truncation_length * truncation_length - zdist * zdist);
2072  if (remdist < 0)
2073  {
2074  continue; // skip if we're too far away
2075  }
2076  }
2077  // sum+=*partial[x][phi][z][ix][iphi][iz] * *q[ix][iphi][iz];
2078  if (r == ir && phi == iphi && z == iz)
2079  {
2080  continue; // dont' compute self-to-self field.
2081  }
2082  sum += Epartial->Get(r - rmin_roi, phi - phimin_roi, z - zmin_roi, ir, iphi, iz) * q->GetChargeInBin(ir, iphi, iz);
2083  }
2084  }
2085  }
2086  // printf("summed field at (%d,%d,%d)=(%f,%f,%f)\n",x,y,z,sum.X(),sum.Y(),sum.Z());
2087  return sum;
2088 }
2089 
2091 {
2092  // do the summation of the inner high-resolution region charges:
2093  //
2094  // bin 0 1 2 ... n-2 n-1
2095  // . . .|.|.|.|.|.|.|. . .
2096  // | | | | | | |
2097  // . . .|.|.|.|.|.|.|. . .
2098  // -----+-+-+-+-+-+-+-----
2099  // . . .|.|.|.|.|.|.|. . .
2100  // -----+-+-+-+-+-+-+-----
2101  // . . .|.|.|.|.|.|.|. . .
2102  // -----+-+-+-+-+-+-+-----
2103  // . . .|.|.|.|.|.|.|. . .
2104  // -----+-+-+-+-+-+-+-----
2105  // . . .|.|.|.|.|.|.|. . .
2106  // | | | | | | |
2107  // . . .|.|.|.|.|.|.|. . .
2108  //
2109  //
2110 
2111  int r_highres_dist = (nr_high - 1) / 2;
2112  int phi_highres_dist = (nphi_high - 1) / 2;
2113  int z_highres_dist = (nz_high - 1) / 2;
2114 
2115  int r_parentlow = floor((r - r_highres_dist) / (r_spacing * 1.0)); // partly enclosed in our high-res region
2116  int r_parenthigh = floor((r + r_highres_dist) / (r_spacing * 1.0)) + 1; // definitely not enclosed in our high-res region
2117  int phi_parentlow = floor((phi - phi_highres_dist) / (phi_spacing * 1.0));
2118  int phi_parenthigh = floor((phi + phi_highres_dist) / (phi_spacing * 1.0)) + 1; // note that this can be bigger than nphi! We keep track of that.
2119  int z_parentlow = floor((z - z_highres_dist) / (z_spacing * 1.0));
2120  int z_parenthigh = floor((z + z_highres_dist) / (z_spacing * 1.0)) + 1;
2121  printf("AnnularFieldSim::sum_local_field_at parents: rlow=%d,philow=%d,zlow=%d,rhigh=%d,phihigh=%d,zhigh=%d\n", r_parentlow, phi_parentlow, z_parentlow, r_parenthigh, phi_parenthigh, z_parenthigh);
2122 
2123  // zero our current qlocal holder:
2124  for (int i = 0; i < q_local->Length(); i++)
2125  {
2126  *(q_local->GetFlat(i)) = 0;
2127  }
2128 
2129  // get the charge involved in the local highres block:
2130  for (int ir = r_parentlow * r_spacing; ir < r_parenthigh * r_spacing; ir++)
2131  {
2132  if (ir < 0)
2133  {
2134  ir = 0;
2135  }
2136  if (ir >= nr)
2137  {
2138  break;
2139  }
2140  int rbin = (ir - r) + r_highres_dist; // index in our highres locale. zeroth bin when we're at max distance below, etc.
2141  if (rbin < 0)
2142  {
2143  rbin = 0;
2144  }
2145  if (rbin >= nr_high)
2146  {
2147  rbin = nr_high - 1;
2148  }
2149  for (int iphi = phi_parentlow * phi_spacing; iphi < phi_parenthigh * phi_spacing; iphi++)
2150  {
2151  // no phi range checks since it's circular.
2152  int phiFilt = FilterPhiIndex(iphi);
2153  int phibin = (iphi - phi) + phi_highres_dist;
2154  if (phibin < 0)
2155  {
2156  phibin = 0;
2157  }
2158  if (phibin >= nphi_high)
2159  {
2160  phibin = nphi_high - 1;
2161  }
2162  for (int iz = z_parentlow * z_spacing; iz < z_parenthigh * z_spacing; iz++)
2163  {
2164  if (iz < 0)
2165  {
2166  iz = 0;
2167  }
2168  if (iz >= nz)
2169  {
2170  break;
2171  }
2172  // if(debugFlag()) printf("%d: AnnularFieldSim::sum_field_at, reading q in f-bin at(r=%d,phi=%d, z=%d)\n",__LINE__,ir,iphi,iz);
2173 
2174  int zbin = (iz - z) + z_highres_dist;
2175  if (zbin < 0)
2176  {
2177  zbin = 0;
2178  }
2179  if (zbin >= nz_high)
2180  {
2181  zbin = nz_high - 1;
2182  }
2183  // printf("filtering in local highres block\n");
2184  q_local->Add(rbin, phibin, zbin, q->GetChargeInBin(ir, phiFilt, iz));
2185  // printf("done filtering in local highres block\n");
2186  }
2187  }
2188  }
2189 
2190  // now q_highres is up to date for our current cell of interest.
2191  // start building our full sum by scaling the local lookup table by q.
2192  // note that the lookup table needs to have already accounted for cell centers.
2193 
2194  TVector3 sum(0, 0, 0);
2195 
2196  // note that Epartial_highres returns zero if we're outside of the global region. q_local will also be zero there.
2197  // these are loops over the position in the epartial highres grid, so relative to the point in question:
2198  // reminder: variables r, phi, and z are global f-bin indices.
2199 
2200  // assuming highres correctly gives a zero when asked for the middle element in each of the last three indices,
2201  // and assuming q_local has no contribution from regions outside the global bounds, this is correct:
2202  for (int ir = 0; ir < nr_high; ir++)
2203  {
2204  for (int iphi = 0; iphi < nphi_high; iphi++)
2205  {
2206  for (int iz = 0; iz < nz_high; iz++)
2207  {
2208  // first three are relative to the roi, last three are relative to the point in the first three. ooph.
2209  if (phi - phimin_roi < 0)
2210  {
2211  printf("%d: Getting with phi=%d\n", __LINE__, phi - phimin_roi);
2212  }
2213  sum += Epartial_highres->Get(r - rmin_roi, phi - phimin_roi, z - zmin_roi, ir, iphi, iz) * q_local->Get(ir, iphi, iz);
2214  }
2215  }
2216  }
2217 
2218  return sum;
2219 }
2220 
2222 {
2223  // now we look for our low_res contribution, which will be the interpolated summed field from the eight blocks closest to our f-bin of interest:
2224 
2225  // find our interpolated position between the eight nearby lowres cells:
2226  // this is similar to the interpolated integral stuff, except we're at integer steps inside integer blocks
2227  // and we have z as well, now.
2228  bool skip[] = {false, false, false, false, false, false, false, false};
2229 
2230  float r0 = r / (1.0 * r_spacing) - 0.5 - rmin_roi_low; // the position in r, in units of r_spacing, starting from the center of the 0th l-bin in the roi.
2231  int r0i = floor(r0); // the integer portion of the position. -- what center is below our position?
2232  float r0d = r0 - r0i; // the decimal portion of the position. -- how far past center are we?
2233  // instead of listing all eight, I'll address these as i%2, (i/2)%2 and (i/4)%2 to avoid typos
2234  int ri[2]; // the r position of the eight cell centers to consider.
2235  ri[0] = r0i;
2236  ri[1] = r0i + 1;
2237  float rw[2]; // the weight of that cell, linear in distance from the center of it
2238  rw[0] = 1 - r0d; // 1 when we're on it, 0 when we're at the other one.
2239  rw[1] = r0d; // 1 when we're on it, 0 when we're at the other one
2240 
2241  // zero out if the requested l-bin is out of the roi (happens if we're close to the outer than the inner edge of the outermost l-bin)
2242  if (ri[0] < 0)
2243  {
2244  for (int i = 0; i < 8; i++)
2245  {
2246  if ((i / 4) % 2 == 0)
2247  {
2248  skip[i] = true; // don't handle contributions from ri[0].
2249  }
2250  }
2251  rw[1] = 1; // and weight like we're dead-center on the outer cells.
2252  }
2253  if (ri[1] >= nr_roi_low)
2254  {
2255  for (int i = 0; i < 8; i++)
2256  {
2257  if ((i / 4) % 2 == 1)
2258  {
2259  skip[i] = true; // don't bother handling ri[1].
2260  }
2261  }
2262  rw[0] = 1; // and weight like we're dead-center on the inner cells.
2263  }
2264 
2265  // now repeat that structure for phi:
2266 
2267  float p0 = phi / (1.0 * phi_spacing) - 0.5 - phimin_roi_low; // the position 'phi' in units of phi_spacing, starting from the center of the 0th bin.
2268  // printf("prepping to filter low, p0=%f, phi=%d, phi_spacing=%d, phimin_roi_low=%d\n",p0,phi,phi_spacing,phimin_roi_low);
2269 
2270  int p0i = floor(p0); // the integer portion of the position. -- what center is below our position?
2271  float p0d = p0 - p0i; // the decimal portion of the position. -- how far past center are we?
2272  int pi[4]; // the local phi position of the four l-bin centers to consider
2273  // printf("filtering low, p0i=%d, nphi_low=%d\n",p0i,nphi_low);
2274  // Q: why do I need to filter here?
2275  // A: Worst case scenario, this produces -1<p0<0. If that wraps around and is still in the roi, I should include it
2276  // A: Likewise, if I get p0>nphi_low, then that means it ought to wrap around and be considered against the other limit.
2277  pi[0] = FilterPhiIndex(p0i, nphi_low);
2278  pi[1] = FilterPhiIndex(p0i + 1, nphi_low);
2279  // having filtered, we will always have a positive number.
2280 
2281  // printf("done filtering low\n");
2282  float pw[2]; // the weight of that cell, linear in distance from the center of it
2283  pw[0] = 1 - p0d; // 1 when we're on it, 0 when we're at the other one.
2284  pw[1] = p0d; // 1 when we're on it, 0 when we're at the other one
2285 
2286  // note that since phi wraps around itself, we have the possibility of being outside by being above/below the opposite end of where we thought we were
2287  // remember these are coordinates wrt the beginning of the roi, and are positive because we filtered. If that rotation didn't put them back in the roi, then they weren't before, either.
2288  if (pi[0] >= nphi_roi_low)
2289  {
2290  for (int i = 0; i < 8; i++)
2291  {
2292  if ((i / 2) % 2 == 0)
2293  {
2294  skip[i] = true; // don't bother handling pi[0].
2295  }
2296  }
2297  pw[1] = 1; // and weight like we're dead-center on the high cells.
2298  }
2299  if (pi[1] >= nphi_roi_low)
2300  {
2301  for (int i = 0; i < 8; i++)
2302  {
2303  if ((i / 2) % 2 == 1)
2304  {
2305  skip[i] = true; // don't bother handling pi[1].
2306  }
2307  }
2308  pw[0] = 1; // and weight like we're dead-center on the high cells.
2309  }
2310 
2311  // and once more for z. ooph.
2312 
2313  float z0 = z / (1.0 * z_spacing) - 0.5 - zmin_roi_low; // the position in r, in units of r_spacing, starting from the center of the 0th bin.
2314  int z0i = floor(z0); // the integer portion of the position. -- what center is below our position?
2315  float z0d = z0 - z0i; // the decimal portion of the position. -- how far past center are we?
2316  // instead of listing all eight, I'll address these as i%2, (i/2)%2 and (i/4)%2 to avoid typos
2317  int zi[2]; // the r position of the eight cell centers to consider.
2318  zi[0] = z0i;
2319  zi[1] = z0i + 1;
2320  float zw[2]; // the weight of that cell, linear in distance from the center of it
2321  zw[0] = 1 - z0d; // 1 when we're on it, 0 when we're at the other one.
2322  zw[1] = z0d; // 1 when we're on it, 0 when we're at the other one
2323 
2324  if (zi[0] < 0)
2325  {
2326  for (int i = 0; i < 8; i++)
2327  {
2328  if ((i) % 2 == 0)
2329  {
2330  skip[i] = true; // don't bother handling zi[0].
2331  }
2332  }
2333  zw[1] = 1; // and weight like we're dead-center on the higher cells.
2334  }
2335  if (zi[1] >= nz_roi_low)
2336  {
2337  for (int i = 0; i < 8; i++)
2338  {
2339  if ((i) % 2 == 1)
2340  {
2341  skip[i] = true; // don't bother handling zi[1].
2342  }
2343  }
2344  zw[0] = 1; // and weight like we're dead-center on the lower cells.
2345  }
2346 
2347  TVector3 sum(0, 0, 0);
2348  // at this point, we should be skipping all destination l-bins that are out-of-bounds.
2349  // note that if any out-of-bounds ones survive somehow, the call to Epartial_lowres will fail loudly.
2350  int lBinEdge[2]; // lower and upper (included) edges of the low-res bin, measured in f-bins, reused per-dimension
2351  int hRegionEdge[2]; // lower and upper edge of the high-res region, measured in f-bins, reused per-dimension.
2352  bool overlapsPhi, overlapsZ; // whether we overlap in R, phi, and z.
2353 
2354  int r_highres_dist = (nr_high - 1) / 2;
2355  int phi_highres_dist = (nphi_high - 1) / 2;
2356  int z_highres_dist = (nz_high - 1) / 2;
2357 
2358  for (int ir = 0; ir < nr_low; ir++)
2359  {
2360  lBinEdge[0] = ir * r_spacing;
2361  lBinEdge[1] = (ir + 1) * r_spacing - 1;
2362  hRegionEdge[0] = r - r_highres_dist;
2363  hRegionEdge[1] = r + r_highres_dist;
2364  bool overlapsR = (lBinEdge[0] <= hRegionEdge[1] && hRegionEdge[0] <= lBinEdge[1]);
2365  for (int iphi = 0; iphi < nphi_low; iphi++)
2366  {
2367  lBinEdge[0] = iphi * phi_spacing;
2368  lBinEdge[1] = (iphi + 1) * phi_spacing - 1;
2369  hRegionEdge[0] = phi - phi_highres_dist;
2370  hRegionEdge[1] = phi + phi_highres_dist;
2371  overlapsPhi = (lBinEdge[0] <= hRegionEdge[1] && hRegionEdge[0] <= lBinEdge[1]);
2372  for (int iz = 0; iz < nz_low; iz++)
2373  {
2374  lBinEdge[0] = iz * z_spacing;
2375  lBinEdge[1] = (iz + 1) * z_spacing - 1;
2376  hRegionEdge[0] = z - z_highres_dist;
2377  hRegionEdge[1] = z + z_highres_dist;
2378  overlapsZ = (lBinEdge[0] <= hRegionEdge[1] && hRegionEdge[0] <= lBinEdge[1]);
2379  // conceptually: see if the l-bin overlaps with the high-res region:
2380  // the high-res region includes all indices from r-(nr_high-1)/2 to r+(nr_high-1)/2.
2381  // each low-res region includes all indices from ir*r_spacing to (ir+1)*r_spacing-1.
2382  if (overlapsR && overlapsPhi && overlapsZ)
2383  {
2384  // if their bounds are interleaved in all dimensions, there is overlap, and we've already summed this region.
2385  continue;
2386  }
2387  else
2388  {
2389  // if(debugFlag()) printf("%d: AnnularFieldSim::sum_field_at, considering l-bin at(r=%d,phi=%d, z=%d)\n",__LINE__,ir,iphi,iz);
2390 
2391  for (int i = 0; i < 8; i++)
2392  {
2393  if (skip[i])
2394  {
2395  continue;
2396  }
2397  if (ri[(i / 4) % 2] + rmin_roi_low == ir && pi[(i / 2) % 2] + phimin_roi_low == iphi && zi[(i) % 2] + zmin_roi_low == iz)
2398  {
2399  printf("considering an l-bins effect on itself, r=%d,phi=%d,z=%d (matches i=%d, not skipped), means we're not interpolating fairly\n", ir, iphi, iz, i);
2400  assert(1 == 2);
2401  }
2402  // the ri, pi, and zi elements are relative to the roi, as needed for Epartial.
2403  // the ir, iphi, and iz are all absolute, as needed for q_lowres
2404  if (pi[(i / 2) % 2] < 0)
2405  {
2406  printf("%d: Getting with phi=%d\n", __LINE__, pi[(i / 2) % 2]);
2407  }
2408  sum += (Epartial_lowres->Get(ri[(i / 4) % 2], pi[(i / 2) % 2], zi[(i) % 2], ir, iphi, iz) * q_lowres->Get(ir, iphi, iz)) * zw[(i) % 2] * pw[(i / 2) % 2] * rw[(i / 4) % 2];
2409  }
2410  }
2411  }
2412  }
2413  }
2414 
2415  return sum;
2416 }
2417 
2419 {
2420  // sum the E field over all nr by ny by nz cells of sources, at the specific position r,phi,z.
2421  // note the specific position in Epartial is in relative coordinates.
2422  // printf("AnnularFieldSim::sum_field_at(r=%d,phi=%d, z=%d)\n",r,phi,z);
2423  TVector3 pos = GetRoiCellCenter(r - rmin_roi, phi - phimin_roi, z - zmin_roi);
2424  TVector3 slicepos = GetRoiCellCenter(r - rmin_roi, 0, z - zmin_roi);
2425  float rotphi = pos.Phi() - slicepos.Phi(); // probably this is phi*step.Phi();
2426 
2427  // caution: if you print progress of each of these, it will print a great deal of data, since this is called per-bin
2428  // unsigned long long totalelements=nr*nphi*nz;
2429  // unsigned long long percent=totalelements/100;
2430 
2431  // unsigned long long el=0;
2432 
2433  TVector3 sum(0, 0, 0);
2434  TVector3 unrotatedField(0, 0, 0);
2435  TVector3 unitField(0, 0, 0);
2436  int phirel;
2437  for (int ir = 0; ir < nr; ir++)
2438  {
2439  for (int iphi = 0; iphi < nphi; iphi++)
2440  {
2441  for (int iz = 0; iz < nz; iz++)
2442  {
2443  // sum+=*partial[x][phi][z][ix][iphi][iz] * *q[ix][iphi][iz];
2444  if (r == ir && phi == iphi && z == iz)
2445  {
2446  continue; // dont' compute self-to-self field.
2447  }
2448  phirel = FilterPhiIndex(iphi - phi);
2449  unitField = Epartial_phislice->Get(r - rmin_roi, 0, z - zmin_roi, ir, phirel, iz);
2450  unitField.RotateZ(rotphi); // previously was rotate by the step.Phi()*phi. //annoying that I can't rename this to 'rotated field' here without unnecessary overhead.
2451 
2452  sum += unitField * q->GetChargeInBin(ir, iphi, iz);
2453  ;
2454 
2455  /*
2456  if(!(el%percent)) {printf("summing phislices %d%%: ",(int)(el/percent));
2457  printf("unit field at (r=%d,p=%d,z=%d) from (ir=%d,ip=%d,iz=%d) is (%E,%E,%E) (xyz), q=%E\n",
2458  r,phi,z,ir,iphi,iz,unitField.X(),unitField.Y(),unitField.Z(),q->GetChargeInBin(ir,iphi,iz));
2459  }
2460  el++;
2461  */
2462  }
2463  }
2464  }
2465  // printf("summed field at (%d,%d,%d)=(%f,%f,%f)\n",x,y,z,sum.X(),sum.Y(),sum.Z());
2466  return sum;
2467 }
2468 
2469 TVector3 AnnularFieldSim::swimToInAnalyticSteps(float zdest, TVector3 start, int steps = 1, int *goodToStep = nullptr)
2470 {
2471  // assume coordinates are given in native units (cm=1 unless that changed!).
2472  double zdist = (zdest - start.Z()) * cm;
2473  double zstep = zdist / steps;
2474  start = start * cm; // scale us to cm.
2475 
2476  TVector3 ret = start;
2477  TVector3 accumulated_distortion(0, 0, 0);
2478  TVector3 accumulated_drift(0, 0, 0);
2479  TVector3 drift_step(0, 0, zstep);
2480 
2481  int rt, pt, zt; // just placeholders for the bounds-checking.
2482  BoundsCase zBound;
2483  for (int i = 0; i < steps; i++)
2484  {
2485  zBound = GetZindexAndCheckBounds(ret.Z(), &zt);
2486  if (zBound == OnLowEdge)
2487  {
2488  // printf("AnnularFieldSIm::swimToInAnalyticSteps requests z-nudge from z=%f to %f\n", ret.Z(), ret.Z()+ALMOST_ZERO);//nudge it in z:
2489  ret.SetZ(ret.Z() + ALMOST_ZERO);
2490  }
2491  if (GetRindexAndCheckBounds(ret.Perp(), &rt) != InBounds || GetPhiIndexAndCheckBounds(FilterPhiPos(ret.Phi()), &pt) != InBounds || (zBound == OutOfBounds))
2492  {
2493  printf(
2494  "AnnularFieldSim::swimToInAnalyticSteps at step %d,"
2495  "asked to swim particle from (%f,%f,%f)(cm) (rphiz)=(%fcm,%frad,%fcm)which is outside the ROI.\n",
2496  i, ret.X() / cm, ret.Y() / cm, ret.Z() / cm, ret.Perp() / cm, ret.Phi(), ret.Z() / cm);
2497  printf(" -- %f <= r < %f \t%f <= phi < %f \t%f <= z < %f \n",
2498  rmin_roi * step.Perp() + rmin, rmax_roi * step.Perp() + rmin,
2499  phimin_roi * step.Phi(), phimax_roi * step.Phi(),
2500  zmin_roi * step.Z(), zmax_roi * step.Z());
2501  printf("Returning last valid position.\n");
2502  if (!(goodToStep == nullptr))
2503  {
2504  *goodToStep = i - 1;
2505  }
2506  return ret * (1.0 / cm);
2507  }
2508  // printf("AnnularFieldSim::swimToInAnalyticSteps at step %d, asked to swim particle from (%f,%f,%f) (rphiz)=(%f,%f,%f).\n",i,ret.X(),ret.Y(),ret.Z(),ret.Perp(),ret.Phi(),ret.Z());
2509  // rcc note: once I put the z distoriton back in, I need to check that ret.Z+zstep is still in bounds:
2510  accumulated_distortion += GetStepDistortion(ret.Z() + zstep, ret, true, true);
2511 
2512  accumulated_drift += drift_step;
2513 
2514  // this seems redundant, but if the distortions are small they may lose precision and stop actually changing the position when step size is small. This allows them to accumulate separately so they can grow properly:
2515  ret = start + accumulated_distortion + accumulated_drift;
2516  }
2517 
2518  return ret * (1.0 / cm);
2519 }
2520 
2521 TVector3 AnnularFieldSim::swimToInSteps(float zdest, const TVector3 &start, int steps = 1, bool interpolate = false, int *goodToStep = nullptr)
2522 {
2523  int *success = nullptr;
2524  TVector3 straightline(start.X(), start.Y(), zdest);
2525  TVector3 distortion = GetTotalDistortion(zdest, start, steps, interpolate, goodToStep, success);
2526  return straightline + distortion;
2527 }
2528 
2529 TVector3 AnnularFieldSim::GetTotalDistortion(float zdest, const TVector3 &start, int steps, bool interpolate, int *goodToStep, int *success)
2530 {
2531  // work in native units is automatic.
2532  // double zdist=(zdest-start.Z())*cm;
2533  // start=start*cm;
2534 
2535  // check the z bounds:
2536  int rt, pt, zt; // just placeholders for the bounds-checking.
2537  BoundsCase zBound;
2538  zBound = GetZindexAndCheckBounds(zdest, &zt);
2539  bool rdrswitch = RdeltaRswitch;
2540  *success = 0;
2541  if (zBound == OutOfBounds)
2542  {
2543  if (hasTwin)
2544  { // check if this destination is valid for our twin, if we have one:
2545  zBound = GetZindexAndCheckBounds(-zdest, &zt);
2546  if (zBound == InBounds || zBound == OnLowEdge)
2547  {
2548  // destination is in the twin
2549  return twin->GetTotalDistortion(zdest, start, steps, interpolate, goodToStep, success);
2550  }
2551  }
2552  // otherwise, we're not in the twin, and default to our usual gripe:
2553  printf("AnnularFieldSim::GetTotalDistortion starting at (%f,%f,%f)=(r%f,p%f,z%f) asked to drift to z=%f, which is outside the ROI. hasTwin= %d. Returning zero_vector.\n", start.X(), start.Y(), start.Z(), start.Perp(), FilterPhiPos(start.Phi()), start.Z(), zdest, (int) hasTwin);
2554  printf(" -- %f <= r < %f \t%f <= phi < %f \t%f <= z < %f \n", rmin_roi * step.Perp() + rmin, rmax_roi * step.Perp() + rmin, phimin_roi * step.Phi(), phimax_roi * step.Phi(), zmin_roi * step.Z(), zmax_roi * step.Z());
2555  *goodToStep = 0;
2556  *success = 0;
2557  return zero_vector; // rcchere
2558  }
2559  else if (zBound == OnLowEdge)
2560  {
2561  // nudge it in z:
2562  zdest += ALMOST_ZERO;
2563  }
2564  zBound = GetZindexAndCheckBounds(start.Z(), &zt);
2565  if (zBound == OutOfBounds)
2566  {
2567  printf("AnnularFieldSim::GetTotalDistortion starting at (%f,%f,%f)=(r%f,p%f,z%f) asked to drift from z=%f, which is outside the ROI. Returning zero_vector.\n", start.X(), start.Y(), start.Z(), start.Perp(), FilterPhiPos(start.Phi()), start.Z(), start.Z());
2568  printf(" -- %f <= r < %f \t%f <= phi < %f \t%f <= z < %f \n", rmin_roi * step.Perp() + rmin, rmax_roi * step.Perp() + rmin, phimin_roi * step.Phi(), phimax_roi * step.Phi(), zmin_roi * step.Z(), zmax_roi * step.Z());
2569  *goodToStep = 0;
2570  *success = 0;
2571  return zero_vector;
2572  }
2573  else if (zBound == OnLowEdge)
2574  {
2575  // nudge it in z:
2576  zdest += ALMOST_ZERO;
2577  }
2578 
2579  // now we are guaranteed the z limits are in range, and don't need to check them again.
2580 
2581  double zstep = (zmax - zmin) / steps; // always a positive number
2582  // printf("Lets do some handy math here, the value of the first zstep: %f \n",zstep);
2583  // printf("Lets do some handy math here, the value of zmax is: %f \n",zmax);
2584  // printf("Lets do some handy math here, the value of zmin is: %f \n",zmin);
2585  // printf("Lets do some handy math here, the value of steps is: %d \n",steps);
2586  if ((zdest - start.Z()) < 0)
2587  {
2588  zstep = -zstep; // It goes two directions, so there will be positive and negative value
2589  }
2590  int integernumbersteps = (zdest - start.Z()) / zstep;
2591  // printf("Lets do some handy math here, the value of integernumbersteps is: %d \n",integernumbersteps);
2592  // printf("Lets do some handy math here, the value of zdest is: %f \n",zdest);
2593  // printf("Lets do some handy math here, the value of start.Z() is: %f \n",start.Z());
2594  // printf("Lets do some handy math here, the value of zstep is: %f \n",zstep);
2595 
2596  // printf("Hello, this is the zmax-zmin length: %f\n", (zmax-zmin));
2597  // printf("Hello, this is the zdest-start.Z(): %f\n", (zdest-start.Z()));
2598  // printf("Hello, this is the zstep number: %f\n", zstep);
2599  // printf("Hello, this is the integernumbersteps number: %d\n", integernumbersteps);
2600 
2601  TVector3 position = start;
2602  TVector3 accumulated_distortion(0, 0, 0);
2603  TVector3 accumulated_drift(0, 0, 0);
2604  TVector3 drift_step(0, 0, zstep);
2605 
2606  // the conceptual approach here is to get the vector distortion in each z step, and use the transverse component of that to update the position of the particle for the next step, while accumulating the total distortion separately from the position. This allows a small residual to accumulate, rather than being lost. We do not correct the z position, so that the stepping does not 'skip' parts of the trajectory.
2607 
2608  // describe what you're doing here:
2609  // printf("Hello again, this is the zdest: %f\n", (zdest));
2610  // printf("Hello again, this is the start.Z(): %f\n", (start.Z()));
2611 
2612  accumulated_distortion += GetStepDistortion(zdest - (integernumbersteps * zstep), position, true, false);
2613 
2614  // short-circuit if there's no travel length:
2615 
2616  position.SetX(start.X() + accumulated_distortion.X());
2617  position.SetY(start.Y() + accumulated_distortion.Y());
2618  position.SetZ(zdest - (integernumbersteps * zstep));
2619 
2620  float StartR, FinalR;
2621  float DeltaR;
2622  for (int i = 0; i < integernumbersteps; i++)
2623  {
2624  if (GetRindexAndCheckBounds(position.Perp(), &rt) != InBounds || GetPhiIndexAndCheckBounds(FilterPhiPos(position.Phi()), &pt) != InBounds || (zBound == OutOfBounds))
2625  {
2626  // printf("AnnularFieldSim::GetTotalDistortion starting at (%f,%f,%f)=(r%f,p%f,z%f) with drift_step=%f, at step %d, asked to swim particle from (%f,%f,%f) (rphiz)=(%f,%f,%f)which is outside the ROI.\n", start.X(), start.Y(), start.Z(), start.Perp(), FilterPhiPos(start.Phi()), start.Z(), zstep, i, position.X(), position.Y(), position.Z(), position.Perp(), position.Phi(), position.Z());
2627  // printf(" -- %f <= r < %f \t%f <= phi < %f \t%f <= z < %f \n", rmin_roi * step.Perp() + rmin, rmax_roi * step.Perp() + rmin, phimin_roi * step.Phi(), phimax_roi * step.Phi(), zmin_roi * step.Z(), zmax_roi * step.Z());
2628  // printf("Returning last good position.\n");
2629  if (!(goodToStep == nullptr))
2630  {
2631  *goodToStep = i - 1, *success = 0;
2632  }
2633  // assert (1==2);
2634  return (accumulated_distortion);
2635  }
2636 
2637  StartR = position.Perp();
2638  // as we accumulate distortions, add these to the x and y positions, but do not change the z position, otherwise we will 'skip' parts of the drift, which is not the intended behavior.
2639  accumulated_distortion += GetStepDistortion(position.Z() + zstep, position, true, false);
2640  position.SetX(start.X() + accumulated_distortion.X());
2641  position.SetY(start.Y() + accumulated_distortion.Y());
2642  FinalR = position.Perp();
2643  // PositionXAfter=start.X() + accumulated_distortion.X();
2644  // PositionYAfter=start.Y() + accumulated_distortion.Y();//
2645 
2646  // StartR=sqrt(PositionXBefore*PositionXBefore+PositionYBefore*PositionYBefore);
2647  DeltaR = FinalR - StartR; // sqrt(PositionXAfter*PositionXAfter+PositionYAfter*PositionYAfter)-StartR;
2648  if (rdrswitch)
2649  {
2650  hRdeltaRComponent->Fill(StartR, DeltaR);
2651  }
2652  position.SetZ(position.Z() + zstep);
2653  }
2654  if (GetRindexAndCheckBounds(position.Perp(), &rt) != InBounds || GetPhiIndexAndCheckBounds(FilterPhiPos(position.Phi()), &pt) != InBounds || (zBound == OutOfBounds))
2655  {
2656  *goodToStep = integernumbersteps - 1;
2657  *success = 0;
2658  return (accumulated_distortion);
2659  }
2660  *goodToStep = integernumbersteps;
2661  *success = 1;
2662  return accumulated_distortion;
2663 }
2664 
2665 void AnnularFieldSim::PlotFieldSlices(const char *filebase, const TVector3 &pos, char which)
2666 {
2667  bool mapEfield = true;
2668  if (which == 'B')
2669  {
2670  mapEfield = false;
2671  }
2672  which = mapEfield ? 'E' : 'B';
2673  char units[5];
2674  if (mapEfield)
2675  {
2676  sprintf(units, "V/cm");
2677  }
2678  else
2679  {
2680  sprintf(units, "T");
2681  }
2682 
2683  printf("plotting field slices for %c field, slicing at (%1.2F,%1.2f,%1.2f)...\n", which, pos.Perp(), FilterPhiPos(pos.Phi()), pos.Z());
2684  std::cout << "file=" << filebase << std::endl;
2685  ;
2686  TString plotfilename = TString::Format("%s.%cfield_slices.pdf", filebase, which);
2687  TVector3 inner = GetInnerEdge();
2688  TVector3 outer = GetOuterEdge();
2689  TVector3 steplocal = GetFieldStep();
2690 
2691  TCanvas *c;
2692 
2693  TH2F *hEfield[3][3];
2694  TH2F *hCharge[3];
2695  TH1F *hEfieldComp[3][3];
2696  char axis[] = "rpzrpz";
2697  float axisval[] = {(float) pos.Perp(), (float) FilterPhiPos(pos.Phi()), (float) pos.Z(), (float) pos.Perp(), (float) FilterPhiPos(pos.Phi()), (float) pos.Z()};
2698  int axn[] = {nr_roi, nphi_roi, nz_roi, nr_roi, nphi_roi, nz_roi};
2699  float axtop[] = {(float) outer.Perp(), 2 * M_PI, (float) outer.Z(), (float) outer.Perp(), 2 * M_PI, (float) outer.Z()};
2700  float axbot[] = {(float) inner.Perp(), 0, (float) inner.Z(), (float) inner.Perp(), 0, (float) inner.Z()};
2701 
2702  // if we are in charge of a twin, extend our axes:
2703  if (hasTwin)
2704  {
2705  // axtop[2]=axtop[5]=(float)(twin->GetOuterEdge().Z());
2706  axbot[2] = axbot[5] = (float) (twin->GetInnerEdge().Z());
2707  }
2708  printf("rpz bounds are %f<r%f\t %f<phi%f\t %f<z%f\n", axbot[0], axtop[0], axbot[1], axtop[1], axbot[2], axtop[2]);
2709  float axstep[6];
2710  for (int i = 0; i < 6; i++)
2711  {
2712  axstep[i] = (axtop[i] - axbot[i]) / (1.0 * axn[i]);
2713  }
2714  TVector3 field;
2715  TVector3 lpos;
2716 
2717  // it's a bit meta, but this loop over axes is a compact way to generate all the histogram titles.
2718  for (int ax = 0; ax < 3; ax++)
2719  {
2720  // loop over which axis slice to take
2721  hCharge[ax] = new TH2F(Form("hCharge%c", axis[ax]),
2722  Form("Spacecharge Distribution in the %c%c plane at %c=%2.3f (C/cm^3);%c;%c",
2723  axis[ax + 1], axis[ax + 2], axis[ax], axisval[ax], axis[ax + 1], axis[ax + 2]),
2724  axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2725  axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2726  for (int i = 0; i < 3; i++)
2727  {
2728  // loop over which axis of the field to read
2729  hEfield[ax][i] = new TH2F(Form("h%cfield%c_%c%c", which, axis[i], axis[ax + 1], axis[ax + 2]),
2730  Form("%c component of %c Field in the %c%c plane at %c=%2.3f (%s);%c;%c",
2731  axis[i], which, axis[ax + 1], axis[ax + 2], axis[ax], axisval[ax], units, axis[ax + 1], axis[ax + 2]),
2732  axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2733  axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2734  hEfieldComp[ax][i] = new TH1F(Form("h%cfieldComp%c_%c%c", which, axis[i], axis[ax + 1], axis[ax + 2]),
2735  Form("Log Magnitude of %c component of %c Field in the %c%c plane at %c=%2.3f;log10(mag)",
2736  axis[i], which, axis[ax + 1], axis[ax + 2], axis[ax], axisval[ax]),
2737  200, -5, 5);
2738  }
2739  }
2740 
2741  float rpz_coord[3];
2742  for (int ax = 0; ax < 3; ax++)
2743  { // we have three sets of 'slices'. the R slice is a 2d plot in phi-z, etc.
2744  rpz_coord[ax] = axisval[ax] + axstep[ax] / 2;
2745  for (int i = 0; i < axn[ax + 1]; i++)
2746  { // for each slice, loop over the bins of the 2d plot:
2747  rpz_coord[(ax + 1) % 3] = axbot[ax + 1] + (i + 0.5) * axstep[ax + 1];
2748  for (int j = 0; j < axn[ax + 2]; j++)
2749  {
2750  rpz_coord[(ax + 2) % 3] = axbot[ax + 2] + (j + 0.5) * axstep[ax + 2];
2751  lpos.SetXYZ(rpz_coord[0], 0, rpz_coord[2]);
2752  lpos.SetPhi(rpz_coord[1]);
2753  if (false && ax == 0)
2754  {
2755  printf("sampling rpz=(%f,%f,%f)=(%f,%f,%f) after conversion to xyz=(%f,%f,%f)\n",
2756  rpz_coord[0], rpz_coord[1], rpz_coord[2],
2757  lpos.Perp(), FilterPhiPos(lpos.Phi()), lpos.Z(), lpos.X(), lpos.Y(), lpos.Z());
2758  }
2759  if (mapEfield)
2760  {
2761  // GetFieldAt automatically asks the twin if we are out of bounds here.
2762  field = GetFieldAt(lpos) * (1.0 * cm / V); // get units so we're drawing in V/cm when we draw.
2763  }
2764  else
2765  {
2766  field = GetBFieldAt(lpos) * (1.0 / Tesla); // get units so we're drawing in Tesla when we draw.
2767  // if (hasTwin && lpos.Z()<0) field=twin->GetBFieldAt(lpos)*1.0/Tesla;
2768  }
2769  field.RotateZ(-rpz_coord[1]); // rotate us so we can read the y component as the phi component
2770  // if (field.Mag()>0) continue; //nothing has mag zero because of the drift field.
2771  hEfield[ax][0]->Fill(rpz_coord[(ax + 1) % 3], rpz_coord[(ax + 2) % 3], field.X());
2772  hEfield[ax][1]->Fill(rpz_coord[(ax + 1) % 3], rpz_coord[(ax + 2) % 3], field.Y());
2773  hEfield[ax][2]->Fill(rpz_coord[(ax + 1) % 3], rpz_coord[(ax + 2) % 3], field.Z());
2774  hCharge[ax]->Fill(rpz_coord[(ax + 1) % 3], rpz_coord[(ax + 2) % 3], GetChargeAt(lpos));
2775  hEfieldComp[ax][0]->Fill((abs(field.X())));
2776  hEfieldComp[ax][1]->Fill((abs(field.Y())));
2777  hEfieldComp[ax][2]->Fill((abs(field.Z())));
2778  }
2779  }
2780  }
2781 
2782  c = new TCanvas("cfieldslices", "electric field", 1200, 800);
2783  c->Divide(4, 3);
2784  gStyle->SetOptStat();
2785  for (int ax = 0; ax < 3; ax++)
2786  {
2787  for (int i = 0; i < 3; i++)
2788  {
2789  c->cd(ax * 4 + i + 1);
2790  gPad->SetRightMargin(0.15);
2791  hEfield[ax][i]->SetStats(false);
2792  hEfield[ax][i]->Draw("colz");
2793  // hEfield[ax][i]->GetListOfFunctions()->Print();
2794  gPad->Modified();
2795  // hEfieldComp[ax][i]->Draw();//"colz");
2796  }
2797  c->cd(ax * 4 + 4);
2798  gPad->SetRightMargin(0.15);
2799  hCharge[ax]->SetStats(false);
2800  hCharge[ax]->Draw("colz");
2801  // pal=dynamic_cast<TPaletteAxis*>(hCharge[ax]->GetListOfFunctions()->FindObject("palette"));
2802  if (false)
2803  {
2804  // pal->SetX1NDC(0.86);
2805  // pal->SetX2NDC(0.91);
2806  gPad->Modified();
2807  }
2808  }
2809  c->SaveAs(plotfilename);
2810  printf("after plotting field slices...\n");
2811  std::cout << "file=" << filebase << std::endl;
2812 
2813  return;
2814 }
2815 
2816 void AnnularFieldSim::GenerateSeparateDistortionMaps(const char *filebase, int nSteps, int r_subsamples, int p_subsamples, int z_subsamples, int /*z_substeps*/, bool andCartesian)
2817 {
2818  // generates the distortion map for one or both sides of the detector, separating them so
2819  // they do not try to interpolate across the CM.
2820 
2821  // 1) pick a map spacing ('s')
2822  TVector3 steplocal(step.Perp() / r_subsamples, 0, step.Z() / z_subsamples);
2823  steplocal.SetPhi(step.Phi() / p_subsamples);
2824  float deltar = steplocal.Perp(); //(rf-ri)/nr;
2825  float deltap = steplocal.Phi(); //(pf-pi)/np;
2826  float deltaz = steplocal.Z(); //(zf-zi)/nz;
2827  bool rdrswitch = RdeltaRswitch;
2828  TVector3 stepzvec(0, 0, deltaz);
2829 
2830  // int nSteps = 500; //how many steps to take in the particle path. hardcoded for now. Think about this later.
2831 
2832  int nSides = 1;
2833  if (hasTwin)
2834  {
2835  nSides = 2;
2836  }
2837  // idea for a faster way to build a map:
2838 
2839  // 2) generate the distortions steplocal.Z() away from the readout
2840  // 3) generate the distortion from (i)*steplocal.Z() away to (i-1) away, then add the interpolated value from the (i-1) layer
2841 
2842  // for interpolation, Henry needs one extra buffer bin on each side.
2843 
2844  // so we define the histogram bounds (the 'h' suffix) to be the full range
2845  // plus an additional step in each direction so interpolation can work at the edges
2846  TVector3 lowerEdge = GetRoiCellCenter(rmin_roi, phimin_roi, zmin_roi);
2847  TVector3 upperEdge = GetRoiCellCenter(rmax_roi - 1, phimax_roi - 1, zmax_roi - 1);
2848  int nph = nphi * p_subsamples + 2; // number of phibins in the histogram
2849  int nrh = nr * r_subsamples + 2; // number of r bins in the histogram
2850  int nzh = nz * z_subsamples + 2; // number of z you get the idea.
2851 
2852  float rih = lowerEdge.Perp() - 0.5 * step.Perp() - steplocal.Perp(); // lower bound of roi, minus one
2853  float rfh = upperEdge.Perp() + 0.5 * step.Perp() + steplocal.Perp(); // upper bound of roi, plus one
2854  float pih = FilterPhiPos(lowerEdge.Phi()) - 0.5 * step.Phi() - steplocal.Phi(); // can't automate this or we'll run afoul of phi looping.
2855  float pfh = FilterPhiPos(upperEdge.Phi()) + 0.5 * step.Phi() + steplocal.Phi(); // can't automate this or we'll run afoul of phi looping.
2856  float zih = lowerEdge.Z() - 0.5 * step.Z() - steplocal.Z(); // lower bound of roi, minus one
2857  float zfh = upperEdge.Z() + 0.5 * step.Z() + steplocal.Z(); // upper bound of roi, plus one
2858  float z_readout = upperEdge.Z() + 0.5 * step.Z(); // readout plane. Note we assume this is positive.
2859 
2860  printf("generating distortion map...\n");
2861  printf("file=%s\n", filebase);
2862  printf("Phi: %d steps from %f to %f (field has %d steps)\n", nph, pih, pfh, GetFieldStepsPhi());
2863  printf("R: %d steps from %f to %f (field has %d steps)\n", nrh, rih, rfh, GetFieldStepsR());
2864  printf("Z: %d steps from %f to %f (field has %d steps)\n", nzh, zih, zfh, GetFieldStepsZ());
2865  TString distortionFilename;
2866  distortionFilename.Form("%s.distortion_map.hist.root", filebase);
2867  TString summaryFilename;
2868  summaryFilename.Form("%s.distortion_summary.pdf", filebase);
2869  TString diffSummaryFilename;
2870  diffSummaryFilename.Form("%s.differential_summary.pdf", filebase);
2871 
2872  TFile *outf = TFile::Open(distortionFilename.Data(), "RECREATE");
2873  outf->cd();
2874 
2875  // actual output maps:
2876 
2877  TH3F *hDistortionR = new TH3F("hDistortionR", "Per-z-bin Distortion in the R direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2878  TH3F *hDistortionP = new TH3F("hDistortionP", "Per-z-bin Distortion in the RPhi direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2879  TH3F *hDistortionZ = new TH3F("hDistortionZ", "Per-z-bin Distortion in the Z direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2880  TH3F *hIntDistortionR = new TH3F("hIntDistortionR", "Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2881  TH3F *hIntDistortionP = new TH3F("hIntDistortionP", "Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2882  TH3F *hIntDistortionZ = new TH3F("hIntDistortionZ", "Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2883 
2884  TH3F *hIntDistortionX = new TH3F("hIntDistortionX", "Integrated X Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2885  TH3F *hIntDistortionY = new TH3F("hIntDistortionY", "Integrated Y Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2886 
2887  const int nMapComponents = 6;
2888  TH3F *hSeparatedMapComponent[2][6]; // side, then xyzrp
2889  TH3F *hValidToStepComponent[2];
2890  TH3F *hReachesReadout[2];
2891  TString side[2];
2892  side[0] = "soloz";
2893  if (hasTwin)
2894  {
2895  side[0] = "posz";
2896  side[1] = "negz";
2897  }
2898  TString sepAxis[] = {"X", "Y", "Z", "R", "P", "RPhi"};
2899  float zlower, zupper;
2900  for (int i = 0; i < nSides; i++)
2901  {
2902  if (i == 0)
2903  { // doing the positive side
2904  zlower = fmin(zih, zfh);
2905  zupper = fmax(zih, zfh);
2906  }
2907  else
2908  {
2909  zlower = -1 * fmax(zih, zfh);
2910  zupper = -1 * fmin(zih, zfh);
2911  }
2912  for (int j = 0; j < nMapComponents; j++)
2913  {
2914  hSeparatedMapComponent[i][j] = new TH3F(Form("hIntDistortion%s_%s", sepAxis[j].Data(), side[i].Data()),
2915  Form("Integrated %s Deflection drifting from (phi,r,z) to z=endcap;phi;r;z (%s side)", sepAxis[j].Data(), side[i].Data()),
2916  nph, pih, pfh, nrh, rih, rfh, nzh, zlower, zupper);
2917  }
2918  hValidToStepComponent[i] = new TH3F(Form("hValidToStep_%s", side[i].Data()),
2919  Form("Step before out of bounds drifting from (phi,r,z) to z=endcap;phi;r;z (%s side)", side[i].Data()),
2920  nph, pih, pfh, nrh, rih, rfh, nzh, zlower, zupper);
2921  hReachesReadout[i] = new TH3F(Form("hReachesReadout_%s", side[i].Data()),
2922  Form("Bool value for checking if successfuly drifting from (phi,r,z) to z=endcap and reaches the readout;phi;r;z (%s side)", side[i].Data()),
2923  nph, pih, pfh, nrh, rih, rfh, nzh, zlower, zupper);
2924  }
2925  if (rdrswitch)
2926  {
2927  hRdeltaRComponent = new TH2F("hRdeltaR", "Bool value of input for a 2D R_versus_deltaR plot from Zstart to Zend;R;deltaR", nrh, rih, rfh, 10 * nrh, rih - 70, rfh + 170);
2928  // hRdeltaRComponent = new TH2F("hRdeltaR","Bool value of input for a 2D R_versus_deltaR plot from Zstart to Zend;R;deltaR", nrh, rih, rfh, 10*nrh, rih-20, rfh-77);
2930  }
2931  // monitor plots, and the position that that plot monitors at:
2932 
2933  // TVector3 pos((nrh/2+0.5)*s.Perp()+rih,0,(nzh/2+0.5)*s.Z()+zih);
2934  TVector3 pos(((int) (nrh / 2) + 0.5) * steplocal.Perp() + rih, 0, zmin + ((int) (nz / 2) + 0.5) * step.Z());
2935  float posphi = ((int) (nph / 2) + 0.5) * steplocal.Phi() + pih;
2936  pos.SetPhi(posphi);
2937  // int xi[3]={nrh/2,nph/2,nzh/2};
2938  int xi[3] = {(int) floor((pos.Perp() - rih) / steplocal.Perp()), (int) floor((posphi - pih) / steplocal.Phi()), (int) floor((pos.Z() - zih) / steplocal.Z())};
2939  if (!hasTwin)
2940  {
2941  printf("rpz slice indices= (%d,%d,%d) (no twin)\n", xi[0], xi[1], xi[2]);
2942  }
2943  int twinz = (-pos.Z() - zih) / steplocal.Z();
2944  if (hasTwin)
2945  {
2946  printf("rpz slice indices= (%d,%d,%d) twinz=%d\n", xi[0], xi[1], xi[2], twinz);
2947  }
2948 
2949  const char axname[] = "rpzrpz";
2950  int axn[] = {nrh, nph, nzh, nrh, nph, nzh};
2951  float axval[] = {(float) pos.Perp(), (float) pos.Phi(), (float) pos.Z(), (float) pos.Perp(), (float) pos.Phi(), (float) pos.Z()};
2952  float axbot[] = {rih, pih, zih, rih, pih, zih};
2953  float axtop[] = {rfh, pfh, zfh, rfh, pfh, zfh};
2954  TH2F *hIntDist[3][3];
2955  TH1F *hRDist[2][3]; // now with a paired friend for the negative side
2956  TH2F *hDiffDist[3][3];
2957  TH1F *hRDiffDist[2][3];
2958  for (int i = 0; i < 3; i++)
2959  {
2960  // loop over which axis of the distortion to read
2961  for (int ax = 0; ax < 3; ax++)
2962  {
2963  // loop over which plane to work in
2964  hDiffDist[ax][i] = new TH2F(TString::Format("hDiffDist%c_%c%c", axname[i], axname[ax + 1], axname[ax + 2]),
2965  TString::Format("%c component of diff. distortion in %c%c plane at %c=%2.3f;%c;%c",
2966  axname[i], axname[ax + 1], axname[ax + 2], axname[ax], axval[ax], axname[ax + 1], axname[ax + 2]),
2967  axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2968  axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2969  hIntDist[ax][i] = new TH2F(TString::Format("hIntDist%c_%c%c", axname[i], axname[ax + 1], axname[ax + 2]),
2970  TString::Format("%c component of int. distortion in %c%c plane at %c=%2.3f;%c;%c",
2971  axname[i], axname[ax + 1], axname[ax + 2], axname[ax], axval[ax], axname[ax + 1], axname[ax + 2]),
2972  axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2973  axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2974  }
2975  hRDist[0][i] = new TH1F(TString::Format("hRDist%c", axname[i]),
2976  TString::Format("%c component of int. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
2977  axname[i], axname[1], axval[1], axname[2], axval[2]),
2978  axn[0], axbot[0], axtop[0]);
2979  hRDist[1][i] = new TH1F(TString::Format("hRDist%cNeg", axname[i]),
2980  TString::Format("%c component of int. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
2981  axname[i], axname[1], axval[1], axname[2], axval[2]),
2982  axn[0], axbot[0], axtop[0]);
2983  hRDiffDist[0][i] = new TH1F(TString::Format("hRDiffDist%c", axname[i]),
2984  TString::Format("%c component of diff. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
2985  axname[i], axname[1], axval[1], axname[2], axval[2]),
2986  axn[0], axbot[0], axtop[0]);
2987  hRDiffDist[1][i] = new TH1F(TString::Format("hRDiffDist%cNeg", axname[i]),
2988  TString::Format("%c component of diff. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
2989  axname[i], axname[1], axval[1], axname[2], axval[2]),
2990  axn[0], axbot[0], axtop[0]);
2991  }
2992 
2993  TVector3 inpart, outpart;
2994  TVector3 diffdistort, distort;
2995  int validToStep;
2996  int successCheck;
2997 
2998  // TTree version:
2999  float partR, partP, partZ;
3000  int ir, ip, iz;
3001  float distortR, distortP, distortZ;
3002  float distortX, distortY;
3003  float diffdistR, diffdistP, diffdistZ;
3004  TTree *dTree = new TTree("dTree", "Distortion per step z");
3005  dTree->Branch("r", &partR);
3006  dTree->Branch("p", &partP);
3007  dTree->Branch("z", &partZ);
3008  dTree->Branch("ir", &ir);
3009  dTree->Branch("ip", &ip);
3010  dTree->Branch("iz", &iz);
3011  dTree->Branch("dr", &distortR);
3012  dTree->Branch("drp", &distortP);
3013  dTree->Branch("dz", &distortZ);
3014 
3015  printf("generating separated distortion map with (%dx%dx%d) grid \n", nrh, nph, nzh);
3016  unsigned long long totalelements = nrh;
3017  totalelements *= nph;
3018  totalelements *= nzh; // breaking up this multiplication prevents a 32bit math overflow
3019  if (hasTwin)
3020  {
3021  totalelements *= 2; // if we have a twin, we have twice as many z bins as we thought.
3022  }
3023 
3024  unsigned long long percent = totalelements / 100;
3025  unsigned long long waypoint = percent * debug_npercent;
3026  printf("total elements = %llu\n", totalelements);
3027 
3028  int el = 0;
3029 
3030  // we want to loop over the entire region to be mapped, but we also need to include
3031  // one additional bin at each edge, to allow the mc drift code to interpolate properly.
3032  // hence we count from -1 to n+1, and manually adjust the position in those edge cases
3033  // to avoid sampling nonphysical regions in r and z. the phi case is free to wrap as
3034  // normal.
3035 
3036  // note that we apply the adjustment to the particle position (inpart) and not the plotted position (partR etc)
3037  inpart.SetXYZ(1, 0, 0);
3038  for (ir = 0; ir < nrh; ir++)
3039  {
3040  partR = (ir + 0.5) * deltar + rih;
3041  if (ir == 0)
3042  {
3043  inpart.SetPerp(partR + deltar);
3044  }
3045  else if (ir == nrh - 1)
3046  {
3047  inpart.SetPerp(partR - deltar);
3048  }
3049  else
3050  {
3051  inpart.SetPerp(partR);
3052  }
3053  for (ip = 0; ip < nph; ip++)
3054  {
3055  partP = (ip + 0.5) * deltap + pih;
3056  inpart.SetPhi(partP);
3057  // since phi loops, there's no need to adjust phis that are out of bounds.
3058  for (iz = 0; iz < nzh; iz++)
3059  {
3060  partZ = (iz) *deltaz + zih; // start us at the EDGE of the bin,
3061  if (iz == 0)
3062  {
3063  inpart.SetZ(partZ + deltaz);
3064  }
3065  else if (iz == nzh - 1)
3066  {
3067  inpart.SetZ(partZ - deltaz);
3068  }
3069  else
3070  {
3071  inpart.SetZ(partZ);
3072  }
3073  partZ += 0.5 * deltaz; // move to center of histogram bin.
3074  for (int localside = 0; localside < nSides; localside++)
3075  {
3076  if (localside == 0)
3077  {
3078  diffdistort = zero_vector; // GetTotalDistortion(inpart.Z() + deltaz, inpart, nSteps, true, &validToStep, &successCheck);
3079  distort = GetTotalDistortion(z_readout, inpart, nSteps, true, &validToStep, &successCheck);
3080  }
3081  else
3082  {
3083  // if we have more than one side,
3084  // flip z coords and do the twin instead:
3085  partZ *= -1; // position to place in histogram
3086  inpart.SetZ(-1 * inpart.Z()); // position to seek in sim
3087  diffdistort = zero_vector; // twin->GetTotalDistortion(inpart.Z() - deltaz, inpart, nSteps, true, &validToStep, &successCheck);
3088  distort = twin->GetTotalDistortion(-z_readout, inpart, nSteps, true, &validToStep, &successCheck);
3089  }
3090 
3091  diffdistort.RotateZ(-inpart.Phi()); // rotate so that distortion components are wrt the x axis
3092  diffdistP = diffdistort.Y(); // the phi component is now the y component.
3093  diffdistR = diffdistort.X(); // and the r component is the x component
3094  diffdistZ = diffdistort.Z();
3095 
3096  distortX = distort.X();
3097  distortY = distort.Y();
3098  distort.RotateZ(-inpart.Phi()); // rotate so that distortion components are wrt the x axis
3099  distortP = distort.Y(); // the phi component is now the y component.
3100  distortR = distort.X(); // and the r component is the x component
3101  distortZ = distort.Z();
3102 
3103  float distComp[nMapComponents]; // by components
3104  distComp[0] = distortX;
3105  distComp[1] = distortY;
3106  distComp[2] = distortZ;
3107  distComp[3] = distortR;
3108  distComp[4] = distortP / partR; // 'P' now refers to phi in radians, rather than unitful
3109  distComp[5] = distortP; // 'RPhi' is the one that correponds to the []meter unitful phi-hat value
3110 
3111  for (int c = 0; c < nMapComponents; c++)
3112  {
3113  hSeparatedMapComponent[localside][c]->Fill(partP, partR, partZ, distComp[c]);
3114  }
3115 
3116  hValidToStepComponent[localside]->Fill(partP, partR, partZ, validToStep);
3117  hReachesReadout[localside]->Fill(partP, partR, partZ, successCheck);
3118 
3119  // recursive integral distortion:
3120  // get others working first!
3121 
3122  // printf("part=(rpz)(%f,%f,%f),distortP=%f\n",partP,partR,partZ,distortP);
3123  hIntDistortionR->Fill(partP, partR, partZ, distortR);
3124  hIntDistortionP->Fill(partP, partR, partZ, distortP);
3125  hIntDistortionZ->Fill(partP, partR, partZ, distortZ);
3126 
3127  if (andCartesian)
3128  {
3129  hIntDistortionX->Fill(partP, partR, partZ, distortX);
3130  hIntDistortionY->Fill(partP, partR, partZ, distortY);
3131  }
3132 
3133  // now we fill particular slices for integral visualizations:
3134  if (ir == xi[0] && localside == 0)
3135  { // r slice
3136  // printf("ir=%d, r=%f (pz)=(%d,%d), distortR=%2.2f, distortP=%2.2f\n",ir,partR,ip,iz,distortR,distortP);
3137  hIntDist[0][0]->Fill(partP, partZ, distortR);
3138  hIntDist[0][1]->Fill(partP, partZ, distortP);
3139  hIntDist[0][2]->Fill(partP, partZ, distortZ);
3140  hDiffDist[0][0]->Fill(partP, partZ, diffdistR);
3141  hDiffDist[0][1]->Fill(partP, partZ, diffdistP);
3142  hDiffDist[0][2]->Fill(partP, partZ, diffdistZ);
3143  }
3144  if (ip == xi[1] && localside == 0)
3145  { // phi slice
3146  // printf("ip=%d, p=%f (rz)=(%d,%d), distortR=%2.2f, distortP=%2.2f\n",ip,partP,ir,iz,distortR,distortP);
3147  hIntDist[1][0]->Fill(partZ, partR, distortR);
3148  hIntDist[1][1]->Fill(partZ, partR, distortP);
3149  hIntDist[1][2]->Fill(partZ, partR, distortZ);
3150  hDiffDist[1][0]->Fill(partZ, partR, diffdistR);
3151  hDiffDist[1][1]->Fill(partZ, partR, diffdistP);
3152  hDiffDist[1][2]->Fill(partZ, partR, diffdistZ);
3153 
3154  if (iz == xi[2] && localside == 0)
3155  { // z slices of phi slices= r line at mid phi, mid z:
3156  hRDist[0][0]->Fill(partR, distortR);
3157  hRDist[0][1]->Fill(partR, distortP);
3158  hRDist[0][2]->Fill(partR, distortZ);
3159  hRDiffDist[0][0]->Fill(partR, diffdistR);
3160  hRDiffDist[0][1]->Fill(partR, diffdistP);
3161  hRDiffDist[0][2]->Fill(partR, diffdistZ);
3162  }
3163  if (hasTwin && iz == twinz && localside == 1)
3164  { // z slices of phi slices= r line at mid phi, mid z:
3165  hRDist[1][0]->Fill(partR, distortR);
3166  hRDist[1][1]->Fill(partR, distortP);
3167  hRDist[1][2]->Fill(partR, distortZ);
3168  hRDiffDist[1][0]->Fill(partR, diffdistR);
3169  hRDiffDist[1][1]->Fill(partR, diffdistP);
3170  hRDiffDist[1][2]->Fill(partR, diffdistZ);
3171  }
3172  }
3173  if (iz == xi[2] && localside == 0)
3174  { // z slice
3175  // printf("iz=%d, z=%f (rp)=(%d,%d), distortR=%2.2f, distortP=%2.2f\n",iz,partZ,ir,ip,distortR,distortP);
3176 
3177  hIntDist[2][0]->Fill(partR, partP, distortR);
3178  hIntDist[2][1]->Fill(partR, partP, distortP);
3179  hIntDist[2][2]->Fill(partR, partP, distortZ);
3180  hDiffDist[2][0]->Fill(partR, partP, diffdistR);
3181  hDiffDist[2][1]->Fill(partR, partP, diffdistP);
3182  hDiffDist[2][2]->Fill(partR, partP, diffdistZ);
3183  }
3184 
3185  if (!(el % waypoint))
3186  {
3187  printf("generating distortions %d%%: ", (int) (el / percent));
3188  printf("distortion at (ir=%d,ip=%d,iz=%d) is (%E,%E,%E)\n",
3189  ir, ip, iz, distortR, distortP, distortZ);
3190  }
3191  el++;
3192  }
3193  }
3194  }
3195  }
3196  printf("Completed distortion generation. Saving outputs...\n");
3197 
3198  TCanvas *canvas = new TCanvas("cdistort", "distortion integrals", 1200, 800);
3199  // take 10 of the bottom of this for data?
3200  printf("was able to make a tcanvas\n");
3201  canvas->cd();
3202  TPad *c = new TPad("cplots", "distortion integral plots", 0, 0.2, 1, 1);
3203  canvas->cd();
3204  TPad *textpad = new TPad("ctext", "distortion integral plots", 0, 0.0, 1, 0.2);
3205  printf("was able to make some tpads\n");
3206 
3207  c->Divide(4, 3);
3208  gStyle->SetOptStat();
3209  printf("was able to interact with gStyle\n");
3210 
3211  for (int i = 0; i < 3; i++)
3212  {
3213  // component
3214  for (int ax = 0; ax < 3; ax++)
3215  {
3216  printf("looping over components i=%d ax=%d\n", i, ax);
3217 
3218  // plane
3219  c->cd(i * 4 + ax + 1);
3220  gPad->SetRightMargin(0.15);
3221  hIntDist[ax][i]->SetStats(false);
3222  hIntDist[ax][i]->Draw("colz");
3223  }
3224 
3225  printf("drawing R profile %d\n", i);
3226 
3227  c->cd(i * 4 + 4);
3228  hRDist[0][i]->SetStats(false);
3229  hRDist[0][i]->SetFillColor(kRed);
3230  hRDist[0][i]->Draw("hist");
3231  if (hasTwin)
3232  {
3233  printf("drawing R profile twin %d\n", i);
3234  hRDist[1][i]->SetStats(false);
3235  hRDist[1][i]->SetLineColor(kBlue);
3236  hRDist[1][i]->Draw("hist,same");
3237  }
3238  }
3239  printf("switching to textpad\n");
3240 
3241  textpad->cd();
3242  float texpos = 0.9;
3243  float texshift = 0.12;
3244  TLatex *tex = new TLatex(0.0, texpos, "Fill Me In");
3245  printf("built TLatex\n");
3246 
3247  tex->SetTextSize(texshift * 0.8);
3248  tex->DrawLatex(0.05, texpos, GetFieldString());
3249  texpos -= texshift;
3250  tex->DrawLatex(0.05, texpos, GetChargeString());
3251  texpos -= texshift;
3252  // tex->DrawLatex(0.05,texpos,Form("Drift Field = %2.2f V/cm",GetNominalE()));texpos-=texshift;
3253  tex->DrawLatex(0.05, texpos, Form("Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3254  texpos -= texshift;
3255  tex->DrawLatex(0.05, texpos, GetLookupString());
3256  texpos -= texshift;
3257  tex->DrawLatex(0.05, texpos, GetGasString());
3258  texpos -= texshift;
3259  if (debug_distortionScale.Mag() > 0)
3260  {
3261  tex->DrawLatex(0.05, texpos, Form("Distortion scaled by (r,p,z)=(%2.2f,%2.2f,%2.2f)", debug_distortionScale.X(), debug_distortionScale.Y(), debug_distortionScale.Z()));
3262  texpos -= texshift;
3263  }
3264  texpos = 0.9;
3265  printf("cd'ing to canvas:\n");
3266 
3267  canvas->cd();
3268  printf("draw1\n");
3269 
3270  c->Draw();
3271  canvas->cd();
3272  printf("draw2\n");
3273  textpad->Draw();
3274  printf("was able to complete drawing on both pads\n");
3275 
3276  canvas->SaveAs(summaryFilename.Data());
3277 
3278  // canvas->cd();
3279  // already done TPad *c=new TPad("cplots","distortion differential plots",0,0.2,1,1);
3280  // canvas->cd();
3281  // already done TPad *textpad=new TPad("ctext","distortion differential plots",0,0.0,1,0.2);
3282  // already done c->Divide(4,3);
3283  // gStyle->SetOptStat();
3284  for (int i = 0; i < 3; i++)
3285  {
3286  // component
3287  for (int ax = 0; ax < 3; ax++)
3288  {
3289  // plane
3290  c->cd(i * 4 + ax + 1);
3291  gPad->SetRightMargin(0.15);
3292  hDiffDist[ax][i]->SetStats(false);
3293  hDiffDist[ax][i]->Draw("colz");
3294  }
3295  c->cd(i * 4 + 4);
3296  hRDiffDist[0][i]->SetStats(false);
3297  hRDiffDist[0][i]->SetFillColor(kRed);
3298  hRDiffDist[0][i]->Draw("hist");
3299  if (hasTwin)
3300  {
3301  hRDiffDist[1][i]->SetStats(false);
3302  hRDiffDist[1][i]->SetLineColor(kBlue);
3303  hRDiffDist[1][i]->Draw("hist same");
3304  }
3305  }
3306  textpad->cd();
3307  texpos = 0.9;
3308  texshift = 0.12;
3309  tex->SetTextSize(texshift * 0.8);
3310  tex->DrawLatex(0.05, texpos, GetFieldString());
3311  texpos -= texshift;
3312  tex->DrawLatex(0.05, texpos, GetChargeString());
3313  texpos -= texshift;
3314  // tex->DrawLatex(0.05,texpos,Form("Drift Field = %2.2f V/cm",GetNominalE()));texpos-=texshift;
3315  tex->DrawLatex(0.05, texpos, Form("Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3316  texpos -= texshift;
3317  tex->DrawLatex(0.05, texpos, GetLookupString());
3318  texpos -= texshift;
3319  tex->DrawLatex(0.05, texpos, GetGasString());
3320  texpos -= texshift;
3321  tex->DrawLatex(0.05, texpos, "Differential Plots");
3322  texpos -= texshift;
3323  if (debug_distortionScale.Mag() > 0)
3324  {
3325  tex->DrawLatex(0.05, texpos, Form("Distortion scaled by (r,p,z)=(%2.2f,%2.2f,%2.2f)", debug_distortionScale.X(), debug_distortionScale.Y(), debug_distortionScale.Z()));
3326  texpos -= texshift;
3327  }
3328  texpos = 0.9;
3329 
3330  canvas->cd();
3331  c->Draw();
3332  canvas->cd();
3333  textpad->Draw();
3334  canvas->SaveAs(diffSummaryFilename.Data());
3335 
3336  printf("saving map histograms to:%s.\n", distortionFilename.Data());
3337 
3338  outf->cd();
3339  for (int i = 0; i < nSides; i++)
3340  {
3341  printf("Saving side '%s'\n", side[i].Data());
3342  for (int j = 0; j < nMapComponents; j++)
3343  {
3344  hSeparatedMapComponent[i][j]->GetSumw2()->Set(0);
3345  hSeparatedMapComponent[i][j]->Write();
3346  }
3347  hValidToStepComponent[i]->GetSumw2()->Set(0);
3348  hValidToStepComponent[i]->Write();
3349 
3350  hReachesReadout[i]->GetSumw2()->Set(0);
3351  hReachesReadout[i]->Write();
3352  }
3353  if (rdrswitch)
3354  {
3355  hRdeltaRComponent->GetSumw2()->Set(0);
3356  hRdeltaRComponent->Write();
3357  }
3358 
3359  hDistortionR->GetSumw2()->Set(0);
3360  hDistortionP->GetSumw2()->Set(0);
3361  hDistortionZ->GetSumw2()->Set(0);
3362  hIntDistortionR->GetSumw2()->Set(0);
3363  hIntDistortionP->GetSumw2()->Set(0);
3364  hIntDistortionZ->GetSumw2()->Set(0);
3365  if (andCartesian)
3366  {
3367  hIntDistortionX->GetSumw2()->Set(0);
3368  hIntDistortionY->GetSumw2()->Set(0);
3369  }
3370  hDistortionR->Write();
3371  hDistortionP->Write();
3372  hDistortionZ->Write();
3373  hIntDistortionR->Write();
3374  hIntDistortionP->Write();
3375  hIntDistortionZ->Write();
3376  if (false && andCartesian)
3377  {
3378  hIntDistortionX->Write();
3379  hIntDistortionY->Write();
3380  }
3381  printf("finished writing histograms\n");
3382  // dTree->Write();
3383  // printf("wrote dTree\n");
3384  outf->Close();
3385  // printf("map:%s.closed\n",distortionFilename.Data());
3386 
3387  printf("wrote separated map and summary to %s.\n", filebase);
3388  printf("map:%s.\n", distortionFilename.Data());
3389  printf("summary:%s.\n", summaryFilename.Data());
3390  // printf("wrote map and summary to %s and %s.\n",distortionFilename.Data(),summaryFilename.Data());
3391  return;
3392 }
3393 
3394 void AnnularFieldSim::GenerateDistortionMaps(const char *filebase, int r_subsamples, int p_subsamples, int z_subsamples, int /*z_substeps*/, bool andCartesian)
3395 {
3396  // generates the distortion map for the full detector instead of one map per side.
3397  // This produces wrong behavior when interpolating across the CM and should not be used
3398  // unless you're doing some debugging/backward compatibility checking. Eventually this should be removed (said in early 2022)
3399  printf("WARNING: You called the version of the distortion generator that generates a unified map. Are you sure you meant to do this?\n");
3400 
3401  // 1) pick a map spacing ('s')
3402  bool makeUnifiedMap = true; // if true, generate a single map across the whole TPC. if false, save two maps, one for each side.
3403 
3404  TVector3 steplocal(step.Perp() / r_subsamples, 0, step.Z() / z_subsamples);
3405  steplocal.SetPhi(step.Phi() / p_subsamples);
3406  float deltar = steplocal.Perp(); //(rf-ri)/nr;
3407  float deltap = steplocal.Phi(); //(pf-pi)/np;
3408  float deltaz = steplocal.Z(); //(zf-zi)/nz;
3409  TVector3 stepzvec(0, 0, deltaz);
3410  int nSteps = 500; // how many steps to take in the particle path. hardcoded for now. Think about this later.
3411 
3412  // idea for a faster way to build a map:
3413 
3414  // 2) generate the distortions steplocal.Z() away from the readout
3415  // 3) generate the distortion from (i)*steplocal.Z() away to (i-1) away, then add the interpolated value from the (i-1) layer
3416 
3417  // for interpolation, Henry needs one extra buffer bin on each side.
3418 
3419  // so we define the histogram bounds (the 'h' suffix) to be the full range
3420  // plus an additional step in each direction so interpolation can work at the edges
3421  TVector3 lowerEdge = GetRoiCellCenter(rmin_roi, phimin_roi, zmin_roi);
3422  TVector3 upperEdge = GetRoiCellCenter(rmax_roi - 1, phimax_roi - 1, zmax_roi - 1);
3423  int nph = nphi * p_subsamples + 2; // nuber of phibins in the histogram
3424  int nrh = nr * r_subsamples + 2; // number of r bins in the histogram
3425  int nzh = nz * z_subsamples + 2; // number of z you get the idea.
3426  int successCheck;
3427 
3428  if (hasTwin && makeUnifiedMap)
3429  { // double the z range if we have a twin. r and phi are the same, unless we had a phi roi...
3430  lowerEdge.SetZ(-1 * upperEdge.Z());
3431  nzh += nz * z_subsamples;
3432  }
3433 
3434  float rih = lowerEdge.Perp() - 0.5 * step.Perp() - steplocal.Perp(); // lower bound of roi, minus one
3435  float rfh = upperEdge.Perp() + 0.5 * step.Perp() + steplocal.Perp(); // upper bound of roi, plus one
3436  float pih = FilterPhiPos(lowerEdge.Phi()) - 0.5 * step.Phi() - steplocal.Phi(); // can't automate this or we'll run afoul of phi looping.
3437  float pfh = FilterPhiPos(upperEdge.Phi()) + 0.5 * step.Phi() + steplocal.Phi(); // can't automate this or we'll run afoul of phi looping.
3438  float zih = lowerEdge.Z() - 0.5 * step.Z() - steplocal.Z(); // lower bound of roi, minus one
3439  float zfh = upperEdge.Z() + 0.5 * step.Z() + steplocal.Z(); // upper bound of roi, plus one
3440  float z_readout = upperEdge.Z() + 0.5 * step.Z(); // readout plane. Note we assume this is positive.
3441 
3442  printf("generating distortion map...\n");
3443  printf("file=%s\n", filebase);
3444  printf("Phi: %d steps from %f to %f (field has %d steps)\n", nph, pih, pfh, GetFieldStepsPhi());
3445  printf("R: %d steps from %f to %f (field has %d steps)\n", nrh, rih, rfh, GetFieldStepsR());
3446  printf("Z: %d steps from %f to %f (field has %d steps)\n", nzh, zih, zfh, GetFieldStepsZ());
3447  TString distortionFilename;
3448  distortionFilename.Form("%s.distortion_map.hist.root", filebase);
3449  TString summaryFilename;
3450  summaryFilename.Form("%s.distortion_summary.pdf", filebase);
3451  TString diffSummaryFilename;
3452  diffSummaryFilename.Form("%s.differential_summary.pdf", filebase);
3453 
3454  TFile *outf = TFile::Open(distortionFilename.Data(), "RECREATE");
3455  outf->cd();
3456 
3457  // actual output maps:
3458 
3459  TH3F *hDistortionR = new TH3F("hDistortionR", "Per-z-bin Distortion in the R direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3460  TH3F *hDistortionP = new TH3F("hDistortionP", "Per-z-bin Distortion in the RPhi direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3461  TH3F *hDistortionZ = new TH3F("hDistortionZ", "Per-z-bin Distortion in the Z direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3462  TH3F *hIntDistortionR = new TH3F("hIntDistortionR", "Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3463  TH3F *hIntDistortionP = new TH3F("hIntDistortionP", "Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3464  TH3F *hIntDistortionZ = new TH3F("hIntDistortionZ", "Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3465 
3466  TH3F *hIntDistortionX = new TH3F("hIntDistortionX", "Integrated X Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3467  TH3F *hIntDistortionY = new TH3F("hIntDistortionY", "Integrated Y Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3468 
3469  /*
3470  TH3F* hNewIntDistortionR=new TH3F("hNewIntDistortionR","Recursively Integrated R Distortion from (r,phi,z) to z=0 (centered in r,phi, and z);phi;r;z",nph,pih,pfh,nrh,rih,rfh,nzh,zih,zfh);
3471  TH3F* hNewIntDistortionP=new TH3F("hNewIntDistortionP","Recursively Integrated R Distortion from (r,phi,z) to z=0 (centered in r,phi, and z);phi;r;z",nph,pih,pfh,nrh,rih,rfh,nzh,zih,zfh);
3472  TH3F* hNewIntDistortionZ=new TH3F("hNewIntDistortionZ","Recursively Integrated R Distortion from (r,phi,z) to z=0 (centered in r,phi, and z);phi;r;z",nph,pih,pfh,nrh,rih,rfh,nzh,zih,zfh);
3473 //for the interchanging distortion maps.
3474  TH2F* hNewIntDistortionR=new TH3F("hNewIntDistortionR","Recursively Integrated R Distortion from (r,phi,z) to z=0 (centered in r,phi, and z);phi;r;z",nph,pih,pfh,nrh,rih,rfh,nzh,zih,zfh);
3475  TH3F* hNewIntDistortionP=new TH3F("hNewIntDistortionP","Recursively Integrated R Distortion from (r,phi,z) to z=0 (centered in r,phi, and z);phi;r;z",nph,pih,pfh,nrh,rih,rfh,nzh,zih,zfh);
3476  TH3F* hNewIntDistortionZ=new TH3F("hNewIntDistortionZ","Recursively Integrated R Distortion from (r,phi,z) to z=0 (centered in r,phi, and z);phi;r;z",nph,pih,pfh,nrh,rih,rfh,nzh,zih,zfh);
3477  */
3478 
3479  // monitor plots, and the position that plot monitors at:
3480 
3481  // TVector3 pos((nrh/2+0.5)*steplocal.Perp()+rih,0,(nzh/2+0.5)*steplocal.Z()+zih);
3482  TVector3 pos(((int) (nrh / 2) + 0.5) * steplocal.Perp() + rih, 0, zmin + ((int) (nz / 2) + 0.5) * step.Z());
3483  float posphi = ((int) (nph / 2) + 0.5) * steplocal.Phi() + pih;
3484  pos.SetPhi(posphi);
3485  // int xi[3]={nrh/2,nph/2,nzh/2};
3486  int xi[3] = {(int) floor((pos.Perp() - rih) / steplocal.Perp()), (int) floor((posphi - pih) / steplocal.Phi()), (int) floor((pos.Z() - zih) / steplocal.Z())};
3487  if (!hasTwin)
3488  {
3489  printf("rpz slice indices= (%d,%d,%d) (no twin)\n", xi[0], xi[1], xi[2]);
3490  }
3491  int twinz = 0; // this is meant to be the matching position to xi[2] in the twin, hence generally -1*pos. Better to just ask the twin rather than trying to calculate it ourselves...
3492  if (hasTwin)
3493  {
3494  twinz = twin->GetZindex(-1 * pos.Z());
3495  printf("rpz slice indices= (%d,%d,%d) twinz=%d\n", xi[0], xi[1], xi[2], twinz);
3496  }
3497 
3498  const char axname[] = "rpzrpz";
3499  int axn[] = {nrh, nph, nzh, nrh, nph, nzh};
3500  float axval[] = {(float) pos.Perp(), (float) pos.Phi(), (float) pos.Z(), (float) pos.Perp(), (float) pos.Phi(), (float) pos.Z()};
3501  float axbot[] = {rih, pih, zih, rih, pih, zih};
3502  float axtop[] = {rfh, pfh, zfh, rfh, pfh, zfh};
3503  TH2F *hIntDist[3][3];
3504  TH1F *hRDist[2][3]; // now with a paired friend for the negative side
3505  TH2F *hDiffDist[3][3];
3506  TH1F *hRDiffDist[2][3];
3507  for (int i = 0; i < 3; i++)
3508  {
3509  // loop over which axis of the distortion to read
3510  for (int ax = 0; ax < 3; ax++)
3511  {
3512  // loop over which plane to work in
3513  hDiffDist[ax][i] = new TH2F(TString::Format("hDiffDist%c_%c%c", axname[i], axname[ax + 1], axname[ax + 2]),
3514  TString::Format("%c component of diff. distortion in %c%c plane at %c=%2.3f;%c;%c",
3515  axname[i], axname[ax + 1], axname[ax + 2], axname[ax], axval[ax], axname[ax + 1], axname[ax + 2]),
3516  axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
3517  axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
3518  hIntDist[ax][i] = new TH2F(TString::Format("hIntDist%c_%c%c", axname[i], axname[ax + 1], axname[ax + 2]),
3519  TString::Format("%c component of int. distortion in %c%c plane at %c=%2.3f;%c;%c",
3520  axname[i], axname[ax + 1], axname[ax + 2], axname[ax], axval[ax], axname[ax + 1], axname[ax + 2]),
3521  axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
3522  axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
3523  }
3524  hRDist[0][i] = new TH1F(TString::Format("hRDist%c", axname[i]),
3525  TString::Format("%c component of int. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
3526  axname[i], axname[1], axval[1], axname[2], axval[2]),
3527  axn[0], axbot[0], axtop[0]);
3528  hRDist[1][i] = new TH1F(TString::Format("hRDist%cNeg", axname[i]),
3529  TString::Format("%c component of int. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
3530  axname[i], axname[1], axval[1], axname[2], axval[2]),
3531  axn[0], axbot[0], axtop[0]);
3532  hRDiffDist[0][i] = new TH1F(TString::Format("hRDist%c", axname[i]),
3533  TString::Format("%c component of diff. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
3534  axname[i], axname[1], axval[1], axname[2], axval[2]),
3535  axn[0], axbot[0], axtop[0]);
3536  hRDiffDist[1][i] = new TH1F(TString::Format("hRDist%cNeg", axname[i]),
3537  TString::Format("%c component of diff. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
3538  axname[i], axname[1], axval[1], axname[2], axval[2]),
3539  axn[0], axbot[0], axtop[0]);
3540  }
3541 
3542  TVector3 inpart, outpart;
3543  TVector3 distort;
3544  int validToStep;
3545 
3546  // TTree version:
3547  float partR, partP, partZ;
3548  int ir, ip, iz;
3549  float distortR, distortP, distortZ;
3550  float distortX, distortY;
3551  float diffdistR, diffdistP, diffdistZ;
3552  TTree *dTree = new TTree("dTree", "Distortion per step z");
3553  dTree->Branch("r", &partR);
3554  dTree->Branch("p", &partP);
3555  dTree->Branch("z", &partZ);
3556  dTree->Branch("ir", &ir);
3557  dTree->Branch("ip", &ip);
3558  dTree->Branch("iz", &iz);
3559  dTree->Branch("dr", &distortR);
3560  dTree->Branch("drp", &distortP);
3561  dTree->Branch("dz", &distortZ);
3562 
3563  printf("generating distortion map with (%dx%dx%d) grid \n", nrh, nph, nzh);
3564  unsigned long long totalelements = nrh;
3565  totalelements *= nph;
3566  totalelements *= nzh; // breaking up this multiplication prevents a 32bit math overflow
3567  unsigned long long percent = totalelements / 100 * debug_npercent;
3568  printf("total elements = %llu\n", totalelements);
3569 
3570  int el = 0;
3571 
3572  // we want to loop over the entire region to be mapped, but we also need to include
3573  // one additional bin at each edge, to allow the mc drift code to interpolate properly.
3574  // hence we count from -1 to n+1, and manually adjust the position in those edge cases
3575  // to avoid sampling nonphysical regions in r and z. the phi case is free to wrap as
3576  // normal.
3577 
3578  // note that we apply the adjustment to the particle position (inpart) and not the plotted position (partR etc)
3579  inpart.SetXYZ(1, 0, 0);
3580  for (ir = 0; ir < nrh; ir++)
3581  {
3582  partR = (ir + 0.5) * deltar + rih;
3583  if (ir == 0)
3584  {
3585  inpart.SetPerp(partR + deltar);
3586  }
3587  else if (ir == nrh - 1)
3588  {
3589  inpart.SetPerp(partR - deltar);
3590  }
3591  else
3592  {
3593  inpart.SetPerp(partR);
3594  }
3595  for (ip = 0; ip < nph; ip++)
3596  {
3597  partP = (ip + 0.5) * deltap + pih;
3598  inpart.SetPhi(partP);
3599  // since phi loops, there's no need to adjust phis that are out of bounds.
3600  for (iz = 0; iz < nzh; iz++)
3601  {
3602  partZ = (iz) *deltaz + zih; // start us at the EDGE of the bin, maybe has problems at the CM when twinned.
3603  if (iz == 0)
3604  {
3605  inpart.SetZ(partZ + deltaz);
3606  }
3607  else if (iz == nzh - 1)
3608  {
3609  inpart.SetZ(partZ - deltaz);
3610  }
3611  else
3612  {
3613  inpart.SetZ(partZ);
3614  }
3615  partZ += 0.5 * deltaz; // move to center of histogram bin.
3616 
3617  // printf("iz=%d, zcoord=%2.2f, bin=%d\n",iz,partZ, hIntDist[0][0]->GetYaxis()->FindBin(partZ));
3618 
3619  // differential distortion:
3620  // be careful with the math of a distortion. The R distortion is NOT the perp() component of outpart-inpart -- that's the transverse magnitude of the distortion!
3621  if (hasTwin && inpart.Z() < 0)
3622  {
3623  distort = twin->GetTotalDistortion(inpart.Z(), inpart + stepzvec, nSteps, true, &validToStep, &successCheck); // step across the cell in the opposite direction, starting at the high side and going to the low side..
3624  }
3625  else
3626  {
3627  distort = GetTotalDistortion(inpart.Z() + deltaz, inpart, nSteps, true, &validToStep, &successCheck);
3628  }
3629  distort.RotateZ(-inpart.Phi()); // rotate so that that is on the x axis
3630  diffdistP = distort.Y(); // the phi component is now the y component.
3631  diffdistR = distort.X(); // and the r component is the x component
3632  diffdistZ = distort.Z();
3633  hDistortionR->Fill(partP, partR, partZ, diffdistR);
3634  hDistortionP->Fill(partP, partR, partZ, diffdistP);
3635  hDistortionZ->Fill(partP, partR, partZ, diffdistZ);
3636  dTree->Fill();
3637 
3638  // integral distortion:
3639  if (hasTwin && makeUnifiedMap && inpart.Z() < 0)
3640  {
3641  distort = twin->GetTotalDistortion(-z_readout, inpart + stepzvec, nSteps, true, &validToStep, &successCheck);
3642  }
3643  else
3644  {
3645  distort = GetTotalDistortion(z_readout, inpart, nSteps, true, &validToStep, &successCheck);
3646  }
3647  distortX = distort.X();
3648  distortY = distort.Y();
3649  distort.RotateZ(-inpart.Phi()); // rotate so that that is on the x axis
3650  distortP = distort.Y(); // the phi component is now the y component.
3651  distortR = distort.X(); // and the r component is the x component
3652  distortZ = distort.Z();
3653 
3654  // recursive integral distortion:
3655  // get others working first!
3656 
3657  // printf("part=(rpz)(%f,%f,%f),distortP=%f\n",partP,partR,partZ,distortP);
3658  hIntDistortionR->Fill(partP, partR, partZ, distortR);
3659  hIntDistortionP->Fill(partP, partR, partZ, distortP);
3660  hIntDistortionZ->Fill(partP, partR, partZ, distortZ);
3661 
3662  if (andCartesian)
3663  {
3664  hIntDistortionX->Fill(partP, partR, partZ, distortX);
3665  hIntDistortionY->Fill(partP, partR, partZ, distortY);
3666  }
3667 
3668  // now we fill particular slices for integral visualizations:
3669  if (ir == xi[0])
3670  { // r slice
3671  // printf("ir=%d, r=%f (pz)=(%d,%d), distortR=%2.2f, distortP=%2.2f\n",ir,partR,ip,iz,distortR,distortP);
3672  hIntDist[0][0]->Fill(partP, partZ, distortR);
3673  hIntDist[0][1]->Fill(partP, partZ, distortP);
3674  hIntDist[0][2]->Fill(partP, partZ, distortZ);
3675  hDiffDist[0][0]->Fill(partP, partZ, diffdistR);
3676  hDiffDist[0][1]->Fill(partP, partZ, diffdistP);
3677  hDiffDist[0][2]->Fill(partP, partZ, diffdistZ);
3678  }
3679  if (ip == xi[1])
3680  { // phi slice
3681  // printf("ip=%d, p=%f (rz)=(%d,%d), distortR=%2.2f, distortP=%2.2f\n",ip,partP,ir,iz,distortR,distortP);
3682  hIntDist[1][0]->Fill(partZ, partR, distortR);
3683  hIntDist[1][1]->Fill(partZ, partR, distortP);
3684  hIntDist[1][2]->Fill(partZ, partR, distortZ);
3685  hDiffDist[1][0]->Fill(partZ, partR, diffdistR);
3686  hDiffDist[1][1]->Fill(partZ, partR, diffdistP);
3687  hDiffDist[1][2]->Fill(partZ, partR, diffdistZ);
3688 
3689  if (iz == xi[2])
3690  { // z slices of phi slices= r line at mid phi, mid z:
3691  hRDist[0][0]->Fill(partR, distortR);
3692  hRDist[0][1]->Fill(partR, distortP);
3693  hRDist[0][2]->Fill(partR, distortZ);
3694  hRDiffDist[0][0]->Fill(partR, diffdistR);
3695  hRDiffDist[0][1]->Fill(partR, diffdistP);
3696  hRDiffDist[0][2]->Fill(partR, diffdistZ);
3697  }
3698  if (hasTwin && iz == twinz)
3699  { // z slices of phi slices= r line at mid phi, mid z:
3700  hRDist[1][0]->Fill(partR, distortR);
3701  hRDist[1][1]->Fill(partR, distortP);
3702  hRDist[1][2]->Fill(partR, distortZ);
3703  hRDiffDist[1][0]->Fill(partR, diffdistR);
3704  hRDiffDist[1][1]->Fill(partR, diffdistP);
3705  hRDiffDist[1][2]->Fill(partR, diffdistZ);
3706  }
3707  }
3708  if (iz == xi[2])
3709  { // z slice
3710  // printf("iz=%d, z=%f (rp)=(%d,%d), distortR=%2.2f, distortP=%2.2f\n",iz,partZ,ir,ip,distortR,distortP);
3711 
3712  hIntDist[2][0]->Fill(partR, partP, distortR);
3713  hIntDist[2][1]->Fill(partR, partP, distortP);
3714  hIntDist[2][2]->Fill(partR, partP, distortZ);
3715  hDiffDist[2][0]->Fill(partR, partP, diffdistR);
3716  hDiffDist[2][1]->Fill(partR, partP, diffdistP);
3717  hDiffDist[2][2]->Fill(partR, partP, diffdistZ);
3718  }
3719 
3720  if (!(el % percent))
3721  {
3722  printf("generating distortions %d%%: ", (int) (debug_npercent * (el / percent)));
3723  printf("distortion at (ir=%d,ip=%d,iz=%d) is (%E,%E,%E)\n",
3724  ir, ip, iz, distortR, distortP, distortZ);
3725  }
3726  el++;
3727  }
3728  }
3729  }
3730 
3731  TCanvas *canvas = new TCanvas("cdistort", "distortion integrals", 1200, 800);
3732  // take 10 of the bottom of this for data?
3733  canvas->cd();
3734  TPad *c = new TPad("cplots", "distortion integral plots", 0, 0.2, 1, 1);
3735  canvas->cd();
3736  TPad *textpad = new TPad("ctext", "distortion integral plots", 0, 0.0, 1, 0.2);
3737  c->Divide(4, 3);
3738  gStyle->SetOptStat();
3739  for (int i = 0; i < 3; i++)
3740  {
3741  // component
3742  for (int ax = 0; ax < 3; ax++)
3743  {
3744  // plane
3745  c->cd(i * 4 + ax + 1);
3746  gPad->SetRightMargin(0.15);
3747  hIntDist[ax][i]->SetStats(false);
3748  hIntDist[ax][i]->Draw("colz");
3749  }
3750  c->cd(i * 4 + 4);
3751  hRDist[0][i]->SetStats(false);
3752  hRDist[0][i]->SetFillColor(kRed);
3753  hRDist[0][i]->Draw("hist");
3754  if (hasTwin)
3755  {
3756  hRDist[1][i]->SetStats(false);
3757  hRDist[1][i]->SetLineColor(kBlue);
3758  hRDist[1][i]->Draw("hist,same");
3759  }
3760  /*
3761  double Cut = 40;
3762  h->SetFillColor(kRed);
3763  TH1F *hNeg = (TH1F*)hRDist[i]->Clone(Form("hNegRDist%d",i));
3764  hNeg->SetFillColor(kGreen);
3765  for (int n = 1; n <= hNeg->GetNbinsX(); n++) {
3766  hNeg->SetBinContent(n,Cut);
3767  }
3768  h3->Draw(); h.Draw("same");
3769  TH1F *h2 = (TH1F*)h->Clone("h2");
3770  h2->SetFillColor(kGray-4);
3771  h2->SetMaximum(Cut);
3772  h2->Draw("same");
3773  */
3774  }
3775  textpad->cd();
3776  float texpos = 0.9;
3777  float texshift = 0.12;
3778  TLatex *tex = new TLatex(0.0, texpos, "Fill Me In");
3779  tex->SetTextSize(texshift * 0.8);
3780  tex->DrawLatex(0.05, texpos, GetFieldString());
3781  texpos -= texshift;
3782  tex->DrawLatex(0.05, texpos, GetChargeString());
3783  texpos -= texshift;
3784  // tex->DrawLatex(0.05,texpos,Form("Drift Field = %2.2f V/cm",GetNominalE()));texpos-=texshift;
3785  tex->DrawLatex(0.05, texpos, Form("Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3786  texpos -= texshift;
3787  tex->DrawLatex(0.05, texpos, GetLookupString());
3788  texpos -= texshift;
3789  tex->DrawLatex(0.05, texpos, GetGasString());
3790  texpos -= texshift;
3791  if (debug_distortionScale.Mag() > 0)
3792  {
3793  tex->DrawLatex(0.05, texpos, Form("Distortion scaled by (r,p,z)=(%2.2f,%2.2f,%2.2f)", debug_distortionScale.X(), debug_distortionScale.Y(), debug_distortionScale.Z()));
3794  texpos -= texshift;
3795  }
3796  texpos = 0.9;
3797 
3798  canvas->cd();
3799  c->Draw();
3800  canvas->cd();
3801  textpad->Draw();
3802  canvas->SaveAs(summaryFilename.Data());
3803 
3804  // canvas->cd();
3805  // already done TPad *c=new TPad("cplots","distortion differential plots",0,0.2,1,1);
3806  // canvas->cd();
3807  // already done TPad *textpad=new TPad("ctext","distortion differential plots",0,0.0,1,0.2);
3808  // already done c->Divide(4,3);
3809  // gStyle->SetOptStat();
3810  for (int i = 0; i < 3; i++)
3811  {
3812  // component
3813  for (int ax = 0; ax < 3; ax++)
3814  {
3815  // plane
3816  c->cd(i * 4 + ax + 1);
3817  gPad->SetRightMargin(0.15);
3818  hDiffDist[ax][i]->SetStats(false);
3819  hDiffDist[ax][i]->Draw("colz");
3820  }
3821  c->cd(i * 4 + 4);
3822  hRDiffDist[0][i]->SetStats(false);
3823  hRDiffDist[0][i]->SetFillColor(kRed);
3824  hRDiffDist[0][i]->Draw("hist");
3825  if (hasTwin)
3826  {
3827  hRDiffDist[1][i]->SetStats(false);
3828  hRDiffDist[1][i]->SetLineColor(kBlue);
3829  hRDiffDist[1][i]->Draw("hist same");
3830  }
3831  }
3832  textpad->cd();
3833  texpos = 0.9;
3834  texshift = 0.12;
3835  tex->SetTextSize(texshift * 0.8);
3836  tex->DrawLatex(0.05, texpos, GetFieldString());
3837  texpos -= texshift;
3838  tex->DrawLatex(0.05, texpos, GetChargeString());
3839  texpos -= texshift;
3840  // tex->DrawLatex(0.05,texpos,Form("Drift Field = %2.2f V/cm",GetNominalE()));texpos-=texshift;
3841  tex->DrawLatex(0.05, texpos, Form("Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3842  texpos -= texshift;
3843  tex->DrawLatex(0.05, texpos, GetLookupString());
3844  texpos -= texshift;
3845  tex->DrawLatex(0.05, texpos, GetGasString());
3846  texpos -= texshift;
3847  tex->DrawLatex(0.05, texpos, "Differential Plots");
3848  texpos -= texshift;
3849  if (debug_distortionScale.Mag() > 0)
3850  {
3851  tex->DrawLatex(0.05, texpos, Form("Distortion scaled by (r,p,z)=(%2.2f,%2.2f,%2.2f)", debug_distortionScale.X(), debug_distortionScale.Y(), debug_distortionScale.Z()));
3852  texpos -= texshift;
3853  }
3854  texpos = 0.9;
3855 
3856  canvas->cd();
3857  c->Draw();
3858  canvas->cd();
3859  textpad->Draw();
3860  canvas->SaveAs(diffSummaryFilename.Data());
3861 
3862  // printf("map:%s.\n",distortionFilename.Data());
3863 
3864  outf->cd();
3865 
3866  hDistortionR->Write();
3867  hDistortionP->Write();
3868  hDistortionZ->Write();
3869  hIntDistortionR->Write();
3870  hIntDistortionP->Write();
3871  hIntDistortionZ->Write();
3872  if (andCartesian)
3873  {
3874  hIntDistortionX->Write();
3875  hIntDistortionY->Write();
3876  }
3877  dTree->Write();
3878  outf->Close();
3879  // printf("map:%s.closed\n",distortionFilename.Data());
3880 
3881  /*
3882  //all histograms associated with this file should be deleted when we closed the file.
3883  //delete the histograms we made:
3884  TH3F *m;
3885  m = (TH3F*)gROOT­>FindObject("hDistortionR"); if(m) { printf("found in memory still.\n"); m­>Delete();}
3886  m = (TH3F*)gROOT­>FindObject("hDistortionP"); if(m) m­>Delete();
3887  m = (TH3F*)gROOT­>FindObject("hDistortionZ"); if(m) m­>Delete();
3888  m = (TH3F*)gROOT­>FindObject("hIntDistortionR"); if(m) m­>Delete();
3889  m = (TH3F*)gROOT­>FindObject("hIntDistortionP"); if(m) m­>Delete();
3890  m = (TH3F*)gROOT­>FindObject("hIntDistortionZ"); if(m) m­>Delete();
3891 
3892  printf("map:%s.deleted via convoluted call. sigh.\n",distortionFilename.Data());
3893  for (int i=0;i<3;i++){
3894  //loop over which axis of the distortion to read
3895  for (int ax=0;ax<3;ax++){
3896  //loop over which plane to work in
3897  hIntDist[ax][i]->Delete();
3898  }
3899  hRDist[i]->Delete();
3900  }
3901  */
3902  printf("wrote map and summary to %s.\n", filebase);
3903  printf("map:%s.\n", distortionFilename.Data());
3904  printf("summary:%s.\n", summaryFilename.Data());
3905  // printf("wrote map and summary to %s and %s.\n",distortionFilename.Data(),summaryFilename.Data());
3906  return;
3907 }
3908 
3909 TVector3 AnnularFieldSim::swimTo(float zdest, const TVector3 &start, bool interpolate, bool useAnalytic)
3910 {
3911  int defaultsteps = 100;
3912  int goodtostep = 0;
3913  if (useAnalytic)
3914  {
3915  return swimToInAnalyticSteps(zdest, start, defaultsteps, &goodtostep);
3916  }
3917  return swimToInSteps(zdest, start, defaultsteps, interpolate, &goodtostep);
3918 }
3919 
3920 TVector3 AnnularFieldSim::GetStepDistortion(float zdest, const TVector3 &start, bool interpolate, bool useAnalytic)
3921 {
3922  // getting the distortion instead of the post-step position allows us to accumulate small deviations from the original position that might be lost in the large number
3923 
3924  // using second order langevin expansion from http://skipper.physics.sunysb.edu/~prakhar/tpc/Papers/ALICE-INT-2010-016.pdf
3925  // TVector3 (*field)[nr][ny][nz]=field_;
3926  int rt, pt, zt; // these are filled by the checkbounds that follow, but are not used.
3927  BoundsCase zBound = GetZindexAndCheckBounds(start.Z(), &zt);
3928  if (GetRindexAndCheckBounds(start.Perp(), &rt) != InBounds || GetPhiIndexAndCheckBounds(FilterPhiPos(start.Phi()), &pt) != InBounds || (zBound != InBounds && zBound != OnHighEdge))
3929  {
3930  // printf("AnnularFieldSim::swimTo asked to swim particle from (xyz)=(%f,%f,%f) which is outside the ROI:\n", start.X(), start.Y(), start.Z());
3931  // printf(" -- %f <= r < %f \t%f <= phi < %f \t%f <= z < %f \n", rmin_roi * step.Perp(), rmax_roi * step.Perp(), phimin_roi * step.Phi(), phimax_roi * step.Phi(), zmin_roi * step.Z(), zmax_roi * step.Z());
3932  // printf("Returning original position.\n");
3933  return start;
3934  }
3935 
3936  // set the direction of the external fields.
3937 
3938  double zdist = zdest - start.Z();
3939 
3940  // short-circuit if there's no travel length:
3941 
3942  if (fabs(zdist) < ALMOST_ZERO * step.Z())
3943  {
3944  // printf("Asked particle from (%f,%f,%f) to z=%f, which is a distance of %fcm. Returning zero.\n", start.X(), start.Y(), start.Z(), zdest, zdist);
3945  return zero_vector;
3946  }
3947 
3948  TVector3 fieldInt; // integral of E field along path
3949  TVector3 fieldIntB; // integral of B field along path
3950 
3951  // note that using analytic takes priority over interpolating todo: clean this up to use a status rther than a pair of flags
3952  if (useAnalytic)
3953  {
3954  fieldInt = analyticFieldIntegral(zdest, start, Efield);
3955  fieldIntB = analyticFieldIntegral(zdest, start, Bfield);
3956  }
3957  else if (interpolate)
3958  {
3959  fieldInt = interpolatedFieldIntegral(zdest, start, Efield);
3960  fieldIntB = interpolatedFieldIntegral(zdest, start, Bfield);
3961  }
3962  else
3963  {
3964  fieldInt = fieldIntegral(zdest, start, Efield);
3965  fieldIntB = fieldIntegral(zdest, start, Bfield);
3966  }
3967 
3968  if (abs(fieldInt.Z() / zdist) < ALMOST_ZERO)
3969  {
3970  printf("GetStepDistortion is attempting to swim with no drift field:\n");
3971  printf("GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
3972  printf("GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
3973  assert(1 == 2);
3974  }
3975  // float fieldz=field_[in3(x,y,0,fx,fy,fz)].Z()+E.Z();// *field[x][y][zi].Z();
3976  double EfieldZ = fieldInt.Z() / zdist; // average field over the path.
3977  double BfieldZ = fieldIntB.Z() / zdist;
3978  // double fieldz=Enominal; // ideal field over path.
3979 
3980  // these values should be with real, not nominal field?
3981  // double mu=abs(vdrift/Enominal);//vdrift in [cm/s], field in [V/cm] hence mu in [cm^2/(V*s)]; should be a positive value. drift velocity over field magnitude, not field direction.
3982  // double omegatau=-mu*Bnominal;//minus sign is for electron charge.
3983  double omegatau = omegatau_nominal; // don't compute this every time!
3984  // or: omegatau=-10*(10*B.Z()/Tesla)*(vdrift/(cm/us))/(fieldz/(V/cm)); //which is the same as my calculation up to a sign.
3985  // printf("omegatau=%f\n",omegatau);
3986 
3987  double T1om = langevin_T1 * omegatau;
3988  double T2om2 = langevin_T2 * omegatau * langevin_T2 * omegatau;
3989  double c0 = 1 / (1 + T2om2); //
3990  double c1 = T1om / (1 + T1om * T1om); // Carlos gets this term wrong. It should be linear on top, quadratic on bottom.
3991  double c2 = T2om2 / (1 + T2om2);
3992 
3993  TVector3 EintOverEz = 1 / EfieldZ * fieldInt; // integral of E/E_z= integral of E / integral of E_z * delta_z
3994  TVector3 BintOverBz = 1 / BfieldZ * fieldIntB; // should this (and the above?) be BfieldZ or Bnominal?
3995 
3996  // really this should be the integral of the ratio, not the ratio of the integrals.
3997  // and should be integrals over the B field, btu for now that's fixed and constant across the region, so not necessary
3998  // there's no reason to do this as r phi. This is an equivalent result, since I handle everything in vectors.
3999  double deltaX = c0 * EintOverEz.X() + c1 * EintOverEz.Y() - c1 * BintOverBz.Y() + c2 * BintOverBz.X();
4000  double deltaY = c0 * EintOverEz.Y() - c1 * EintOverEz.X() + c1 * BintOverBz.X() + c2 * BintOverBz.Y();
4001  // strictly, for deltaZ we want to integrate v'(E)*(E-E0)dz and v''(E)*(E-E0)^2 dz, but over a short step the field is constant, and hence this can be a product of the integral and not an integral of the product:
4002 
4003  double vprime = (5000 * cm / s) / (100 * V / cm); // hard-coded value for 50:50. Eventually this needs to be part of the constructor, as do most of the repeated math terms here.
4004  double vdoubleprime = 0; // neglecting the v'' term for now. Fair? It's pretty linear at our operating point, and it would require adding an additional term to the field integral.
4005 
4006  // note: as long as my step is very small, I am essentially just reading the field at a point and multiplying by the step size.
4007  // hence integral of P dz, where P is a function of the fields, int|P(E(x,y,z))dz=P(int|E(x,y,z)dz/deltaZ)*deltaZ
4008  // hence: , eg, int|E^2dz=(int|Edz)^2/deltaz
4009  double deltaZ = vprime / vdrift * (fieldInt.Z() - zdist * Enominal) + vdoubleprime / vdrift * (fieldInt.Z() - Enominal * zdist) * (fieldInt.Z() - Enominal * zdist) / (2 * zdist) - 0.5 / zdist * (EintOverEz.X() * EintOverEz.X() + EintOverEz.Y() * EintOverEz.Y()) + c1 / zdist * (EintOverEz.X() * BintOverBz.Y() - EintOverEz.Y() * BintOverBz.X()) + c2 / zdist * (EintOverEz.X() * BintOverBz.X() + EintOverEz.Y() * BintOverBz.Y()) + c2 / zdist * (BintOverBz.X() * BintOverBz.X() + BintOverBz.Y() * BintOverBz.Y()); // missing v'' term.
4010 
4011  if (false)
4012  {
4013  printf("GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4014  printf("GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4015  printf("GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4016  printf("GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
4017  printf("GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4018  printf("GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4019  }
4020 
4021  // if (abs(deltaZ / zdist) > 0.25)
4022  if (false)
4023  {
4024  printf("GetStepDistortion produced a very large zdistortion!\n");
4025  printf("GetStepDistortion: zdist=%2.4f, deltaZ=%2.4f, Delta/z=%2.4f\n", zdist, deltaZ, deltaZ / zdist);
4026  printf("GetStepDistortion: average field Z: Ez=%2.4fV/cm, Bz=%2.4fT\n", EfieldZ / (V / cm), BfieldZ / Tesla);
4027  printf("GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4028  printf("GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4029  printf("GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4030  printf("GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
4031  printf("GetStepDistortion: EfieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4032  printf("GetStepDistortion: BfieldInt=(%E,%E,%E)\n", fieldIntB.X(), fieldIntB.Y(), fieldIntB.Z());
4033  printf("GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4034  // assert(false);
4035  }
4036 
4037  if (abs(deltaX) < 1E-20 && !(chargeCase == NoSpacecharge))
4038  {
4039  printf("GetStepDistortion produced a very small deltaX: %E\n", deltaX);
4040  printf("GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4041  printf("GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4042  printf("GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4043  printf("GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
4044  printf("GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4045  printf("GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4046  // assert(1==2);
4047  }
4048 
4049  // if (!(abs(deltaX) < 1E3))
4050  if (false)
4051  {
4052  printf("GetStepDistortion produced a very large deltaX: %E\n", deltaX);
4053  printf("GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4054  printf("GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4055  printf("GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4056  printf("GetStepDistortion: (%2.4f,%2.4f,%2.4f) (rp)=(%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), start.Perp(), start.Phi(), zdest);
4057  printf("GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4058  printf("GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4059  assert(1 == 2);
4060  }
4061 
4062  // deltaZ=0;//temporary removal.
4063 
4064  TVector3 shift(deltaX, deltaY, deltaZ);
4065  if (debug_distortionScale.Mag() > 0)
4066  { // debug code to scale the resulting distortions
4067  shift.RotateZ(-start.Phi());
4068  // TVector3 localScale=debug_distortionScale;
4069  // localScale.RotateZ(start.Phi());
4070  shift.SetXYZ(shift.X() * debug_distortionScale.X(), shift.Y() * debug_distortionScale.Y(), shift.Z() * debug_distortionScale.Z());
4071  shift.RotateZ(start.Phi());
4072  }
4073 
4074  return shift;
4075 }
4076 
4077 // putting all the getters here out of the way:
4079 {
4080  if (lookupCase == LookupCase::Full3D)
4081  {
4082  return Form("Full3D (%d x %d x %d) with (%d x %d x %d) roi", nr, nphi, nz, nr_roi, nphi_roi, nz_roi);
4083  }
4084 
4085  if (lookupCase == LookupCase::PhiSlice)
4086  {
4087  return Form("PhiSlice (%d x %d x %d) with (%d x 1 x %d) roi", nr, nphi, nz, nr_roi, nz_roi);
4088  }
4089 
4090  return "broken";
4091 }
4093 {
4094  return Form("vdrift=%2.2fcm/us, Enom=%2.2fV/cm, Bnom=%2.2fT, omtau=%2.4E", vdrift / (cm / us), Enominal / (V / cm), Bnominal / Tesla, omegatau_nominal);
4095 }
4096 
4098 {
4099  return Form("%s, %s", Efieldname.c_str(), Bfieldname.c_str());
4100 }
4101 
4102 TVector3 AnnularFieldSim::GetFieldAt(const TVector3 &pos)
4103 {
4104  // assume pos is in native units (see header)
4105 
4106  int r, p, z;
4107 
4108  if (GetRindexAndCheckBounds(pos.Perp(), &r) == BoundsCase::OutOfBounds)
4109  {
4110  return zero_vector;
4111  }
4112  if (GetPhiIndexAndCheckBounds(FilterPhiPos(pos.Phi()), &p) == BoundsCase::OutOfBounds)
4113  {
4114  return zero_vector;
4115  }
4116  if (GetZindexAndCheckBounds(pos.Z(), &z) == BoundsCase::OutOfBounds)
4117  {
4118  if (hasTwin)
4119  {
4120  return twin->GetFieldAt(pos);
4121  }
4122  return zero_vector;
4123  }
4124  return Efield->Get(r, p, z);
4125 }
4126 
4127 TVector3 AnnularFieldSim::GetBFieldAt(const TVector3 &pos)
4128 {
4129  // assume pos is in native units (see header)
4130 
4131  int r, p, z;
4132 
4133  if (GetRindexAndCheckBounds(pos.Perp(), &r) == BoundsCase::OutOfBounds)
4134  {
4135  return zero_vector;
4136  }
4137  if (GetPhiIndexAndCheckBounds(FilterPhiPos(pos.Phi()), &p) == BoundsCase::OutOfBounds)
4138  {
4139  return zero_vector;
4140  }
4141  if (GetZindexAndCheckBounds(pos.Z(), &z) == BoundsCase::OutOfBounds)
4142  {
4143  if (hasTwin)
4144  {
4145  return twin->GetBFieldAt(pos);
4146  }
4147  return zero_vector;
4148  }
4149  return Bfield->Get(r, p, z);
4150 }
4151 
4152 float AnnularFieldSim::GetChargeAt(const TVector3 &pos)
4153 {
4154  int z;
4155  BoundsCase zbound = GetZindexAndCheckBounds(pos.Z(), &z); //==BoundsCase::OutOfBounds) return zero_vector;
4156  if (zbound == OutOfBounds)
4157  {
4158  if (hasTwin)
4159  {
4160  return twin->GetChargeAt(pos);
4161  }
4162  printf("Caution: tried to read charge at zbin=%d! No twin available to handle this\n", z);
4163  return -999;
4164  }
4165 
4166  return q->GetChargeAtPosition(pos.Perp(), FilterPhiPos(pos.Phi()), pos.Z()); // because tvectors take position to be -phi to phi, we always have to filter.
4167  // actually, we should probably just yield to that assumption in more places to speed this up.
4168  /*
4169  //assume pos is in native units (see header)
4170  int r, p, z;
4171 
4172  //get the bounds, but we don't want to actually check the cases, because the charge can go outside the vector region.
4173  r = GetRindex(pos.Perp());
4174  p = GetPhiIndex(pos.Phi());
4175 
4176  BoundsCase zbound = GetZindexAndCheckBounds(pos.Z(), &z); //==BoundsCase::OutOfBounds) return zero_vector;
4177  if (zbound == OutOfBounds && hasTwin) return twin->GetChargeAt(pos);
4178  return q->Get(r, p, z);
4179  */
4180 }