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