17 #include <boost/format.hpp>
22 #define ALMOST_ZERO 0.00001
25 int r,
int roi_r0,
int roi_r1,
int ,
int ,
26 int phi,
int roi_phi0,
int roi_phi1,
int ,
int ,
27 int z,
int roi_z0,
int roi_z1,
int ,
int ,
31 if (roi_r0 >= r || roi_r1 > r || roi_r0 >= roi_r1)
35 if (roi_phi0 >= phi || roi_phi1 > phi || roi_phi0 >= roi_phi1)
37 printf(
"phi roi is out of range or spans the wrap-around. Please spare me that math.\n");
40 if (roi_z0 >= z || roi_z1 > z || roi_z0 >= roi_z1)
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);
72 rmin = in_innerRadius *
cm;
73 rmax = in_outerRadius *
cm;
99 printf(
"AnnularFieldSim::AnnularFieldSim set variables nr=%d, nphi=%d, nz=%d\n",
nr,
nphi,
nz);
105 step.SetXYZ(1, 0, 0);
113 for (
int i = 0; i <
q->
Length(); i++)
124 printf(
"AnnularFieldSim::AnnularFieldSim set roi variables rmin=%d phimin=%d zmin=%d rmax=%d phimax=%d zmax=%d\n",
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.0
e6));
177 printf(
"lookupCase==HybridRes\n");
187 printf(
"lookupCase==PhiSlice\n");
210 printf(
"lookupCase==Analytic (or NoLookup)\n");
232 printf(
"Ran into wrong lookupCase logic in constructor.\n");
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,
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,
249 printf(
"AnnularFieldSim::OldConstructor building AnnularFieldSim with ChargeCase::FromFile\n");
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,
258 r, roi_r0, roi_r1, 1, 3,
259 phi, roi_phi0, roi_phi1, 1, 3,
260 z, roi_z0, roi_z1, 1, 3,
265 printf(
"AnnularFieldSim::OldConstructor building AnnularFieldSim with local_size=1 in all dimensions, lowres_spacing=1 in all dimensions\n");
276 printf(
"AnnularFieldSim::OldConstructor building AnnularFieldSim with roi=full in all dimensions\n");
286 printf(
"something's asking for 'at' with r=%f, which is index=%d\n", at.Perp(), r_position);
294 TVector3 field(0, 0, 0);
301 TVector3 delr = at - from;
310 field.SetMag(
k_perm * 1 / (delr * delr));
318 double delphi =
abs(at.DeltaPhi(from));
322 double Er =
green->
Er(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
327 Ephi =
green->
Ephi(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
329 double Ez =
green->
Ez(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
330 field.SetXYZ(-Er, -Ephi, -Ez);
332 field.RotateZ(at.Phi());
353 if (range < 0) range =
nphi;
365 if (p >= range || p < 0)
367 printf(
"AnnularFieldSim::FilterPhiIndex asked to filter %d, which is more than range=%d out of bounds. Check what called this.\n", phi, range);
375 float r0f = (pos -
rmin) /
step.Perp();
382 float p0f = (
pos) /
step.Phi();
383 int phitemp = floor(p0f);
390 float z0f = (pos -
zmin) /
step.Z();
399 float r0f = (pos -
rmin) /
step.Perp();
425 float p0f = (
pos) /
step.Phi();
426 int phitemp = floor(p0f);
459 float z0f = (pos -
zmin) /
step.Z();
500 if (!rOkay || !phiOkay)
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);
506 int dir = (start.Z() > zdest ? -1 : 1);
527 bool startOkay = (startBound ==
InBounds);
530 if (!startOkay || !endOkay)
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);
541 return dir * integral;
559 if (!rOkay || !phiOkay)
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);
565 int dir = (start.Z() > zdest ? -1 : 1);
586 bool startOkay = (startBound ==
InBounds);
589 if (!startOkay || !endOkay)
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);
595 TVector3 fieldInt(0, 0, 0);
598 for (
int i = zi; i < zf; i++)
602 fieldInt += tf *
step.Z();
616 return dir * fieldInt;
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);
649 float ravg = (r0 +
r1) / 2.0 + 0.5;
650 if (phi0 > phi1) phi1 +=
nphi;
653 printf(
"phi1(%d)<=phi0(%d) even after boosting phi1. check what called this!\n", phi1, phi0);
656 float phiavg = (r0 +
r1) / 2.0 + 0.5;
659 float zavg = (z0 +
z1) / 2.0 + 0.5;
662 c.SetPerp((ravg) *
step.Perp() +
rmin);
663 c.SetPhi((phiavg) *
step.Phi());
676 float rout = rin +
step.Perp();
678 float rMid = (4 * sin(
step.Phi() / 2) * (pow(rout, 3) - pow(rin, 3)) / (3 *
step.Phi() * (pow(rout, 2) - pow(rin, 2))));
680 c.SetPhi((phi + 0.5) *
step.Phi());
681 c.SetZ((z + 0.5) *
step.Z());
690 float r0 = (start.Perp() -
rmin) /
step.Perp() - 0.5;
692 float r0d = r0 - r0i;
695 ri[2] = ri[3] = r0i + 1;
697 rw[0] = rw[1] = 1 - r0d;
700 bool skip[] = {
false,
false,
false,
false};
715 float p0 = (start.Phi()) /
step.Phi() - 0.5;
717 float p0d = p0 - p0i;
722 pw[0] = pw[2] = 1 - p0d;
739 int dir = (start.Z() < zdest ? 1 : -1);
763 if (!startOkay || !endOkay)
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);
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++)
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]);
788 TVector3 fieldInt(0, 0, 0), partialInt;
790 for (
int i = 0; i < 4; i++)
797 partialInt.SetXYZ(0, 0, 0);
798 for (
int j = zi; j < zf; j++)
815 fieldInt += rw[i] * pw[i] * partialInt;
818 return dir * fieldInt;
830 double ifc_radius = 83.5 *
cm;
831 double ofc_radius = 254.5 *
cm;
832 double tpc_halfz = 250 *
cm;
835 double totalcharge = 0;
836 double localcharge = 0;
839 for (
int ifr = 0; ifr <
nr; ifr++)
841 for (
int ifphi = 0; ifphi <
nphi; ifphi++)
843 for (
int ifz = 0; ifz <
nz; ifz++)
846 pos = pos * (1.0 /
cm);
850 totalcharge += localcharge;
851 q->
Add(ifr, ifphi, ifz, localcharge);
855 printf(
"AnnularFieldSim::load_analytic_spacecharge: Total charge Q=%E Coulombs\n", totalcharge);
865 for (
int ifr = 0; ifr <
nr; ifr++)
868 for (
int ifphi = 0; ifphi <
nphi; ifphi++)
871 for (
int ifz = 0; ifz <
nz; ifz++)
886 TFile fieldFile(filename.c_str(),
"READ");
887 TTree *fTree = (TTree *) (fieldFile.Get(treename.c_str()));
888 Efieldname =
"E:" + filename +
":" + treename;
892 fTree->SetBranchAddress(
"r", &r);
893 fTree->SetBranchAddress(
"er", &fr);
894 fTree->SetBranchAddress(
"z", &z);
895 fTree->SetBranchAddress(
"ez", &fz);
908 TFile fieldFile(filename.c_str(),
"READ");
909 TTree *fTree = (TTree *) (fieldFile.Get(treename.c_str()));
910 Bfieldname =
"B:" + filename +
":" + treename;
913 fTree->SetBranchAddress(
"r", &r);
914 fTree->SetBranchAddress(
"br", &fr);
915 fTree->SetBranchAddress(
"z", &z);
916 fTree->SetBranchAddress(
"bz", &fz);
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);
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);
952 bool phiSymmetry = (phiptr == 0);
953 int lowres_factor = 10;
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);
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);
960 for (
int i = 0; i < 3; i++)
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);
966 int nEntries = source->GetEntries();
967 for (
int i = 0; i < nEntries; i++)
970 float zval = *zptr * zsign;
975 htEntries->Fill(*phiptr, *rptr, zval);
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);
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);
986 for (
int j = 0; j <
nphi; j++)
988 htEntries->Fill(j *
step.Phi(), *rptr, zval);
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);
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);
1001 for (
int i = 0; i <
nphi; i++)
1003 for (
int j = 0; j <
nr; j++)
1005 for (
int k = 0;
k <
nz;
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)
1018 (*field)->Set(j, i,
k, fieldvec);
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++)
1027 for (
int j = 0; j <
nr; j++)
1029 for (
int k = 0;
k <
nz;
k++)
1032 int bin = htEntries->FindBin(
FilterPhiPos(cellcenter.Phi()), cellcenter.Perp(), cellcenter.Z());
1033 if (htEntries->GetBinContent(bin) == 0)
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)
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());
1046 (*field)->Set(j, i,
k, fieldvec);
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;
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;
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();
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);
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;
1115 if (!isChargeDensity)
1117 for (
int i = 0; i < hphin; i++)
1119 float phi = hphimin + hphistep * i;
1120 for (
int j = 0; j < hrn; j++)
1122 float r = hrmin + hrstep * j;
1123 for (
int k = 0;
k < hzn;
k++)
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);
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;
1140 for (
int i = 0; i < new_nphi; i++)
1142 float phi = hphimin + new_phistep * (i + 0.5);
1143 float hphi = (phi - hphimin) / hphistep;
1144 bool phisafe = ((hphi - 0.75) > 0) && ((hphi + 0.75) < hphin);
1145 for (
int j = 0; j < new_nr; j++)
1147 float r = hrmin + new_rstep * (j + 0.5);
1148 float hr = (r - hrmin) / hrstep;
1149 bool rsafe = ((hr - 0.5) > 0) && ((hr + 0.5) < hrn);
1150 for (
int k = 0;
k < new_nz;
k++)
1152 float z = hzmin + new_zstep * (
k + 0.5);
1153 float hz = (z - hzmin) / hzstep;
1154 bool zsafe = ((hz - 0.5) > 0) && ((hz + 0.5) < hzn);
1156 if (phisafe && rsafe && zsafe)
1159 resampled->Fill(phi, r, z, hist->Interpolate(phi, r, z));
1163 int bin = hist->FindBin(phi, r, z);
1164 resampled->Fill(phi, r, z, hist->GetBinContent(bin));
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;
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);
1198 float hrstep = (hrmax - hrmin) / hrn;
1199 float hphistep = (hphimax - hphimin) / hphin;
1200 float hzstep = (hzmax - hzmin) / hzn;
1203 int hnzmin = (
zmin - hzmin) / hzstep;
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);
1208 for (
int i = 0; i <
q->
Length(); i++)
1219 double totalcharge = 0;
1222 int localr, localphi, localz;
1224 for (
int i = (
rmin - hrmin) / hrstep; i < hrn; i++)
1226 hr = hrmin + hrstep * (i + 0.5);
1227 localr = (hr -
rmin) /
step.Perp();
1231 printf(
"Loading from histogram has r out of range! r=%f < rmin=%f\n", hr,
rmin);
1236 printf(
"Loading from histogram has r out of range! r=%f > rmax=%f\n", hr,
rmax);
1240 for (
int j = 0; j < hphin; j++)
1242 hphi = hphimin + hphistep * (j + 0.5);
1243 localphi = hphi /
step.Phi();
1244 if (localphi >=
nphi)
1253 for (
int k = hnzmin;
k < hnzmax;
k++)
1255 hz = hzmin + hzstep * (
k + 0.5);
1256 localz = (hz + zoffset -
zmin) /
step.Z();
1274 if (isChargeDensity)
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;
1281 qbin = chargescale * hist->GetBinContent(hist->GetBin(j + 1, i + 1,
k + 1)) *
C;
1284 totalcharge += qbin;
1286 q->
Add(localr, localphi, localz, qbin);
1291 printf(
"AnnularFieldSim::load_spacecharge: Total charge Q=%E Coulombs\n", totalcharge /
C);
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);
1304 for (
int ifr = 0; ifr <
nr; ifr++)
1307 for (
int ifphi = 0; ifphi <
nphi; ifphi++)
1310 for (
int ifz = 0; ifz <
nz; ifz++)
1324 int rcell, phicell, zcell;
1331 if (rb == BoundsCase::OutOfBounds || pb == BoundsCase::OutOfBounds || zb == BoundsCase::OutOfBounds)
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);
1335 if (rcell < 0 || phicell < 0 || zcell < 0)
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);
1341 q->
Add(rcell, phicell, zcell, coulombs *
C);
1344 printf(
"write the code you didn't write earlier about reloading the lowres map.\n");
1383 unsigned long long totalelements =
nr_roi;
1386 unsigned long long percent = totalelements / 100 *
debug_npercent;
1387 printf(
"total elements = %llu\n", totalelements *
nr *
nphi *
nz);
1399 if (!(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());
1425 printf(
"lookupCase==Full3D\n");
1431 printf(
"lookupCase==HybridRes\n");
1437 printf(
"Populating lookup: lookupCase==PhiSlice\n");
1442 printf(
"Populating lookup: lookupCase==Analytic ===> skipping!\n");
1446 printf(
"Populating lookup: lookupCase==NoLookup ===> skipping!\n");
1460 printf(
"populating full lookup table for (%dx%dx%d)x(%dx%dx%d) grid\n",
1465 totalelements *=
nr;
1466 totalelements *=
nphi;
1467 totalelements *=
nz;
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);
1482 for (
int ior = 0; ior <
nr; ior++)
1484 for (
int iophi = 0; iophi <
nphi; iophi++)
1486 for (
int ioz = 0; ioz <
nz; ioz++)
1489 if (!(el % percent))
printf(
"populate_full3d_lookup %d%%\n", (
int) (
debug_npercent * el / percent));
1494 if (ifr == ior && ifphi == iophi && ifz == ioz)
1513 TVector3
at(1, 0, 0);
1514 TVector3 from(1, 0, 0);
1515 TVector3
zero(0, 0, 0);
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;
1523 static int nfbinsin[3][3][3];
1524 for (
int i = 0; i < 3; i++)
1526 for (
int j = 0; j < 3; j++)
1528 for (
int k = 0;
k < 3;
k++)
1530 nfbinsin[i][j][
k] = 0;
1538 TVector3 currentf, newf;
1541 int r_parentlow = floor((ifr - r_highres_dist) / (
r_spacing * 1.0));
1542 int r_parenthigh = floor((ifr + r_highres_dist) / (
r_spacing * 1.0)) + 1;
1543 int r_startpoint = r_parentlow *
r_spacing;
1544 int r_endpoint = r_parenthigh *
r_spacing;
1549 bool phi_parentlow_wrapped = (ifphi - phi_highres_dist < 0);
1551 if (phi_parentlow_wrapped) phi_startpoint -=
nphi;
1553 int phi_parenthigh = floor(
FilterPhiIndex(ifphi + phi_highres_dist) / (phi_spacing * 1.0)) + 1;
1554 bool phi_parenthigh_wrapped = (ifphi + phi_highres_dist >=
nphi);
1556 if (phi_parenthigh_wrapped) phi_endpoint +=
nphi;
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;
1565 int z_endpoint = z_parenthigh *
z_spacing;
1576 for (
int ir = r_startpoint; ir < r_endpoint; ir++)
1581 if (ir >=
nr)
break;
1583 int rbin = (ir - ifr) + r_highres_dist;
1596 for (
int iphi = phi_startpoint; iphi < phi_endpoint; iphi++)
1600 int phibin = (iphi - ifphi) + phi_highres_dist;
1612 for (
int iz = z_startpoint; iz < z_endpoint; iz++)
1615 if (iz >=
nz)
break;
1616 int zbin = (iz - ifz) + z_highres_dist;
1631 nfbinsin[rcell][phicell][zcell]++;
1632 int nf = nfbinsin[rcell][phicell][zcell];
1637 if (zcell != 1 || rcell != 1 || phicell != 1)
1641 if (iphi_rel < 0)
printf(
"%d: Getting with phi=%d\n", __LINE__, iphi_rel);
1645 newf = (currentf * (nf - 1) +
calc_unit_field(at, from)) * (1 / (nf * 1.0));
1652 if (ifr == rbin && ifphi == phibin && ifz == zbin)
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;
1677 int r_low, r_high, phi_low, phi_high, z_low, z_high;
1683 fr_high = fr_low + r_spacing - 1;
1684 if (fr_high >=
nr) fr_high =
nr - 1;
1688 fphi_high = fphi_low + phi_spacing - 1;
1689 if (fphi_high >=
nphi) fphi_high =
nphi - 1;
1693 fz_high = fz_low + z_spacing - 1;
1694 if (fz_high >=
nz) fz_high =
nz - 1;
1699 for (
int ior = 0; ior <
nr_low; ior++)
1702 r_high = r_low + r_spacing - 1;
1705 if (r_high >=
nr) r_high =
nr - 1;
1706 for (
int iophi = 0; iophi <
nphi_low; iophi++)
1709 phi_high = phi_low + phi_spacing - 1;
1710 if (phi_high >=
nphi) phi_high =
nphi - 1;
1712 for (
int ioz = 0; ioz <
nz_low; ioz++)
1715 z_high = z_low + z_spacing - 1;
1716 if (z_high >=
nz) z_high =
nz - 1;
1720 if (ifr == ior && ifphi == iophi && ifz == ioz)
1747 unsigned long long totalelements =
nr;
1748 totalelements *=
nphi;
1749 totalelements *=
nz;
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);
1764 for (
int ior = 0; ior <
nr; ior++)
1766 for (
int iophi = 0; iophi <
nphi; iophi++)
1768 for (
int ioz = 0; ioz <
nz; ioz++)
1774 if (ifr == ior && 0 == iophi && ifz == ioz)
1776 if (!(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());
1789 if (!(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());
1810 unsigned long long totalelements =
nr;
1811 totalelements *=
nphi;
1812 totalelements *=
nz;
1815 unsigned long long percent = totalelements / 100 *
debug_npercent;
1816 printf(
"total elements = %llu\n", totalelements);
1818 TFile *input = TFile::Open(sourcefile,
"READ");
1820 TTree *tInfo = (TTree *) (input->Get(
"info"));
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);
1839 printf(
"param\tobj\tfile\n");
1840 printf(
"nr\t%d\t%d\n",
nr, file_nr);
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);
1852 if (file_rmin !=
rmin || file_rmax !=
rmax ||
1853 file_zmin !=
zmin || file_zmax !=
zmax ||
1856 file_nr !=
nr || file_np !=
nphi || file_nz !=
nz)
1858 printf(
"file parameters do not match fieldsim parameters:\n");
1863 TTree *tLookup = (TTree *) (input->Get(
"phislice"));
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);
1871 tLookup->SetBranchAddress(
"iz_source", &ioz);
1872 tLookup->SetBranchAddress(
"iz_target", &ifz);
1873 tLookup->SetBranchAddress(
"Evec", &unitf);
1876 printf(
"%s has %lld entries\n", sourcefile, tLookup->GetEntries());
1877 for (
int i = 0; i < (
int) totalelements; i++)
1880 tLookup->GetEntry(i);
1884 if (!(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());
1899 unsigned long long totalelements =
nr;
1900 totalelements *=
nphi;
1901 totalelements *=
nz;
1904 unsigned long long percent = totalelements / 100 *
debug_npercent;
1905 printf(
"total elements = %llu\n", totalelements);
1907 TFile *output = TFile::Open(destfile,
"RECREATE");
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);
1923 printf(
"info tree built.\n");
1925 TTree *tLookup =
new TTree(
"phislice",
"Phislice Lookup Table");
1926 int ior, ifr, iophi, ioz, ifz;
1928 tLookup->Branch(
"ir_source", &ior);
1929 tLookup->Branch(
"ir_target", &ifr);
1930 tLookup->Branch(
"ip_source", &iophi);
1932 tLookup->Branch(
"iz_source", &ioz);
1933 tLookup->Branch(
"iz_target", &ifz);
1934 tLookup->Branch(
"Evec", &unitf);
1935 printf(
"lookup tree built.\n");
1942 for (ior = 0; ior <
nr; ior++)
1944 for (iophi = 0; iophi <
nphi; iophi++)
1946 for (ioz = 0; ioz <
nz; ioz++)
1952 if (!(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());
1977 printf(
"AnnularFieldSim::setFlatFields(B=%f Tesla,E=%f V/cm)\n", B, E);
1980 sprintf(fieldstr,
"%f", E);
1981 Efieldname =
"E:Flat:" + std::string(fieldstr);
1982 sprintf(fieldstr,
"%f", B);
1983 Bfieldname =
"B:Flat:" + std::string(fieldstr);
2009 TVector3
sum(0, 0, 0);
2032 if (
debugFlag())
printf(
"summed field at (%d,%d,%d)=(%f,%f,%f)\n", r, phi, z, sum.X(), sum.Y(), sum.Z());
2042 TVector3
sum(0, 0, 0);
2043 float rdist, phidist, zdist, remdist;
2044 for (
int ir = 0; ir <
nr; ir++)
2048 rdist =
abs(ir - r);
2050 if (remdist < 0)
continue;
2052 for (
int iphi = 0; iphi <
nphi; iphi++)
2056 phidist = fmin(
abs(iphi - phi),
abs(
abs(iphi - phi) - nphi));
2058 if (remdist < 0)
continue;
2060 for (
int iz = 0; iz <
nz; iz++)
2064 zdist =
abs(iz - z);
2066 if (remdist < 0)
continue;
2069 if (r == ir && phi == iphi && z == iz)
continue;
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;
2103 int r_parentlow = floor((r - r_highres_dist) / (
r_spacing * 1.0));
2104 int r_parenthigh = floor((r + r_highres_dist) / (
r_spacing * 1.0)) + 1;
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;
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);
2119 if (ir >=
nr)
break;
2120 int rbin = (ir -
r) + r_highres_dist;
2121 if (rbin < 0) rbin = 0;
2127 int phibin = (iphi -
phi) + phi_highres_dist;
2128 if (phibin < 0) phibin = 0;
2133 if (iz >=
nz)
break;
2136 int zbin = (iz -
z) + z_highres_dist;
2137 if (zbin < 0) zbin = 0;
2150 TVector3
sum(0, 0, 0);
2158 for (
int ir = 0; ir <
nr_high; ir++)
2160 for (
int iphi = 0; iphi <
nphi_high; iphi++)
2162 for (
int iz = 0; iz <
nz_high; iz++)
2181 bool skip[] = {
false,
false,
false,
false,
false,
false,
false,
false};
2184 int r0i = floor(r0);
2185 float r0d = r0 - r0i;
2197 for (
int i = 0; i < 8; i++)
2198 if ((i / 4) % 2 == 0)
2204 for (
int i = 0; i < 8; i++)
2205 if ((i / 4) % 2 == 1)
2215 int p0i = floor(p0);
2216 float p0d = p0 - p0i;
2235 for (
int i = 0; i < 8; i++)
2236 if ((i / 2) % 2 == 0)
2242 for (
int i = 0; i < 8; i++)
2243 if ((i / 2) % 2 == 1)
2251 int z0i = floor(z0);
2252 float z0d = z0 - z0i;
2263 for (
int i = 0; i < 8; i++)
2270 for (
int i = 0; i < 8; i++)
2276 TVector3
sum(0, 0, 0);
2281 bool overlapsR, overlapsPhi, overlapsZ;
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;
2287 for (
int ir = 0; ir <
nr_low; ir++)
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++)
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++)
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]);
2311 if (overlapsR && overlapsPhi && overlapsZ)
2320 for (
int i = 0; i < 8; i++)
2322 if (skip[i])
continue;
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);
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];
2348 float rotphi = pos.Phi() - slicepos.Phi();
2358 TVector3
sum(0, 0, 0);
2359 TVector3 unrotatedField(0, 0, 0);
2360 TVector3 unitField(0, 0, 0);
2362 for (
int ir = 0; ir <
nr; ir++)
2364 for (
int iphi = 0; iphi <
nphi; iphi++)
2366 for (
int iz = 0; iz <
nz; iz++)
2369 if (r == ir && phi == iphi && z == iz)
continue;
2372 unitField.RotateZ(rotphi);
2374 sum += unitField *
q->
Get(ir, iphi, iz);
2394 double zdist = (zdest - start.Z()) *
cm;
2395 double zstep = zdist /
steps;
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);
2405 for (
int i = 0; i <
steps; i++)
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",
2423 printf(
"Returning last valid position.\n");
2424 if (!(goodToStep == 0)) *goodToStep = i - 1;
2425 return ret * (1.0 /
cm);
2430 accumulated_drift += drift_step;
2433 ret = start + accumulated_distortion + accumulated_drift;
2436 return ret * (1.0 /
cm);
2441 TVector3 straightline(start.X(), start.Y(), zdest);
2443 return straightline + distortion;
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);
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());
2492 double zstep = (zdest - start.Z()) / steps;
2495 TVector3 accumulated_distortion(0, 0, 0);
2496 TVector3 accumulated_drift(0, 0, 0);
2497 TVector3 drift_step(0, 0, zstep);
2501 for (
int i = 0; i <
steps; i++)
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());
2508 printf(
"Returning last good position.\n");
2509 if (!(goodToStep == 0)) *goodToStep = i - 1;
2511 return (accumulated_distortion);
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);
2520 *goodToStep =
steps;
2521 return accumulated_distortion;
2526 bool mapEfield =
true;
2531 which = mapEfield ?
'E' :
'B';
2535 sprintf(units,
"V/cm");
2539 sprintf(units,
"T");
2542 printf(
"plotting field slices for %c field...\n", which);
2543 std::cout <<
"file=" << filebase << std::endl;
2545 TString plotfilename = TString::Format(
"%s.%cfield_slices.pdf", filebase.c_str(), which);
2552 TH2F *hEfield[3][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()};
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()};
2569 for (
int i = 0; i < 6; i++)
2571 axstep[i] = (axtop[i] - axbot[i]) / (1.0 * axn[i]);
2577 for (
int ax = 0;
ax < 3;
ax++)
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++)
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]),
2601 for (
int ax = 0;
ax < 3;
ax++)
2603 rpz_coord[
ax] = axisval[
ax] + axstep[
ax] / 2;
2604 for (
int i = 0; i < axn[
ax + 1]; i++)
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++)
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]);
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());
2627 field.RotateZ(-rpz_coord[1]);
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())));
2640 c =
new TCanvas(
"cfieldslices",
"electric field", 1200, 800);
2642 gStyle->SetOptStat();
2643 for (
int ax = 0;
ax < 3;
ax++)
2645 for (
int i = 0; i < 3; i++)
2647 c->cd(
ax * 4 + i + 1);
2648 gPad->SetRightMargin(0.15);
2649 hEfield[
ax][i]->SetStats(0);
2650 hEfield[
ax][i]->Draw(
"colz");
2656 gPad->SetRightMargin(0.15);
2657 hCharge[
ax]->SetStats(0);
2658 hCharge[
ax]->Draw(
"colz");
2667 c->SaveAs(plotfilename);
2668 printf(
"after plotting field slices...\n");
2669 std::cout <<
"file=" << filebase << std::endl;
2677 TVector3
s(
step.Perp() / r_subsamples, 0,
step.Z() / z_subsamples);
2678 s.SetPhi(
step.Phi() / p_subsamples);
2679 float deltar =
s.Perp();
2680 float deltap =
s.Phi();
2681 float deltaz =
s.Z();
2682 TVector3 stepzvec(0, 0, deltaz);
2698 int nph =
nphi * p_subsamples + 2;
2699 int nrh =
nr * r_subsamples + 2;
2700 int nzh =
nz * z_subsamples + 2;
2702 float rih = lowerEdge.Perp() - 0.5 *
step.Perp() -
s.Perp();
2703 float rfh = upperEdge.Perp() + 0.5 *
step.Perp() +
s.Perp();
2706 float zih = lowerEdge.Z() - 0.5 *
step.Z() -
s.Z();
2707 float zfh = upperEdge.Z() + 0.5 *
step.Z() +
s.Z();
2708 float z_readout = upperEdge.Z() + 0.5 *
step.Z();
2710 printf(
"generating distortion map...\n");
2711 printf(
"file=%s\n", filebase);
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);
2722 TFile *outf = TFile::Open(distortionFilename.Data(),
"RECREATE");
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);
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);
2737 TH3F *hSeparatedMapComponent[2][5];
2745 char sepAxis[] =
"XYZRP";
2746 float zlower, zupper;
2747 for (
int i = 0; i < nSides; i++)
2751 zlower = fmin(zih, zfh);
2752 zupper = fmax(zih, zfh);
2756 zlower = -1 * fmax(zih, zfh);
2757 zupper = -1 * fmin(zih, zfh);
2759 for (
int j = 0; j < 5; j++)
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);
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;
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);
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];
2786 TH2F *hDiffDist[3][3];
2787 TH1F *hRDiffDist[2][3];
2788 for (
int i = 0; i < 3; i++)
2791 for (
int ax = 0;
ax < 3;
ax++)
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]);
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]);
2823 TVector3 inpart, outpart;
2824 TVector3 diffdistort, distort;
2828 float partR, partP, partZ;
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);
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;
2848 unsigned long long percent = totalelements / 100 *
debug_npercent;
2849 printf(
"total elements = %llu\n", totalelements);
2860 inpart.SetXYZ(1, 0, 0);
2861 for (ir = 0; ir < nrh; ir++)
2863 partR = (ir + 0.5) * deltar + rih;
2866 inpart.SetPerp(partR + deltar);
2868 else if (ir == nrh - 1)
2870 inpart.SetPerp(partR - deltar);
2874 inpart.SetPerp(partR);
2876 for (ip = 0; ip < nph; ip++)
2878 partP = (ip + 0.5) * deltap + pih;
2879 inpart.SetPhi(partP);
2881 for (iz = 0; iz < nzh; iz++)
2883 partZ = (iz) *deltaz + zih;
2886 inpart.SetZ(partZ + deltaz);
2888 else if (iz == nzh - 1)
2890 inpart.SetZ(partZ - deltaz);
2896 partZ += 0.5 * deltaz;
2897 for (
int side = 0; side < nSides; side++)
2901 diffdistort =
GetTotalDistortion(inpart.Z() + deltaz, inpart, nSteps,
true, &validToStep);
2909 inpart.SetZ(-1 * inpart.Z());
2914 diffdistort.RotateZ(-inpart.Phi());
2915 diffdistP = diffdistort.Y();
2916 diffdistR = diffdistort.X();
2917 diffdistZ = diffdistort.Z();
2919 distortX = distort.X();
2920 distortY = distort.Y();
2921 distort.RotateZ(-inpart.Phi());
2922 distortP = distort.Y();
2923 distortR = distort.X();
2924 distortZ = distort.Z();
2927 distComp[0] = distortX;
2928 distComp[1] = distortY;
2929 distComp[2] = distortZ;
2930 distComp[3] = distortR;
2931 distComp[4] = distortP;
2933 for (
int c = 0;
c < 5;
c++)
2935 hSeparatedMapComponent[side][
c]->Fill(partP, partR, partZ, distComp[
c]);
2942 hIntDistortionR->Fill(partP, partR, partZ, distortR);
2943 hIntDistortionP->Fill(partP, partR, partZ, distortP);
2944 hIntDistortionZ->Fill(partP, partR, partZ, distortZ);
2948 hIntDistortionX->Fill(partP, partR, partZ, distortX);
2949 hIntDistortionY->Fill(partP, partR, partZ, distortY);
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);
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);
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);
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);
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);
3004 if (!(el % percent))
3007 printf(
"distortion at (ir=%d,ip=%d,iz=%d) is (%E,%E,%E)\n",
3008 ir, ip, iz, distortR, distortP, distortZ);
3016 TCanvas *canvas =
new TCanvas(
"cdistort",
"distortion integrals", 1200, 800);
3019 TPad *
c =
new TPad(
"cplots",
"distortion integral plots", 0, 0.2, 1, 1);
3021 TPad *textpad =
new TPad(
"ctext",
"distortion integral plots", 0, 0.0, 1, 0.2);
3023 gStyle->SetOptStat();
3024 for (
int i = 0; i < 3; i++)
3027 for (
int ax = 0;
ax < 3;
ax++)
3030 c->cd(i * 4 +
ax + 1);
3031 gPad->SetRightMargin(0.15);
3032 hIntDist[
ax][i]->SetStats(0);
3033 hIntDist[
ax][i]->Draw(
"colz");
3036 hRDist[0][i]->SetStats(0);
3037 hRDist[0][i]->SetFillColor(kRed);
3038 hRDist[0][i]->Draw(
"hist");
3041 hRDist[1][i]->SetStats(0);
3042 hRDist[1][i]->SetLineColor(kBlue);
3043 hRDist[1][i]->Draw(
"hist,same");
3048 float texshift = 0.12;
3049 TLatex *tex =
new TLatex(0.0, texpos,
"Fill Me In");
3050 tex->SetTextSize(texshift * 0.8);
3056 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3073 canvas->SaveAs(summaryFilename.Data());
3081 for (
int i = 0; i < 3; i++)
3084 for (
int ax = 0;
ax < 3;
ax++)
3087 c->cd(i * 4 +
ax + 1);
3088 gPad->SetRightMargin(0.15);
3089 hDiffDist[
ax][i]->SetStats(0);
3090 hDiffDist[
ax][i]->Draw(
"colz");
3093 hRDiffDist[0][i]->SetStats(0);
3094 hRDiffDist[0][i]->SetFillColor(kRed);
3095 hRDiffDist[0][i]->Draw(
"hist");
3098 hRDiffDist[1][i]->SetStats(0);
3099 hRDiffDist[1][i]->SetLineColor(kBlue);
3100 hRDiffDist[1][i]->Draw(
"hist same");
3106 tex->SetTextSize(texshift * 0.8);
3112 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3118 tex->DrawLatex(0.05, texpos,
"Differential Plots");
3131 canvas->SaveAs(diffSummaryFilename.Data());
3136 for (
int i = 0; i < nSides; i++)
3138 for (
int j = 0; j < 5; j++)
3140 hSeparatedMapComponent[i][j]->GetSumw2()->Set(0);
3141 hSeparatedMapComponent[i][j]->Write();
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);
3152 hIntDistortionX->GetSumw2()->Set(0);
3153 hIntDistortionY->GetSumw2()->Set(0);
3155 hDistortionR->Write();
3156 hDistortionP->Write();
3157 hDistortionZ->Write();
3158 hIntDistortionR->Write();
3159 hIntDistortionP->Write();
3160 hIntDistortionZ->Write();
3163 hIntDistortionX->Write();
3164 hIntDistortionY->Write();
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());
3180 bool makeUnifiedMap =
true;
3182 TVector3
s(
step.Perp() / r_subsamples, 0,
step.Z() / z_subsamples);
3183 s.SetPhi(
step.Phi() / p_subsamples);
3184 float deltar =
s.Perp();
3185 float deltap =
s.Phi();
3186 float deltaz =
s.Z();
3187 TVector3 stepzvec(0, 0, deltaz);
3201 int nph =
nphi * p_subsamples + 2;
3202 int nrh =
nr * r_subsamples + 2;
3203 int nzh =
nz * z_subsamples + 2;
3205 if (
hasTwin && makeUnifiedMap)
3207 lowerEdge.SetZ(-1 * upperEdge.Z());
3208 nzh +=
nz * z_subsamples;
3211 float rih = lowerEdge.Perp() - 0.5 *
step.Perp() -
s.Perp();
3212 float rfh = upperEdge.Perp() + 0.5 *
step.Perp() +
s.Perp();
3215 float zih = lowerEdge.Z() - 0.5 *
step.Z() -
s.Z();
3216 float zfh = upperEdge.Z() + 0.5 *
step.Z() +
s.Z();
3217 float z_readout = upperEdge.Z() + 0.5 *
step.Z();
3219 printf(
"generating distortion map...\n");
3220 printf(
"file=%s\n", filebase);
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);
3231 TFile *outf = TFile::Open(distortionFilename.Data(),
"RECREATE");
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);
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);
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;
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);
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];
3275 TH2F *hDiffDist[3][3];
3276 TH1F *hRDiffDist[2][3];
3277 for (
int i = 0; i < 3; i++)
3280 for (
int ax = 0;
ax < 3;
ax++)
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]);
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]);
3312 TVector3 inpart, outpart;
3317 float partR, partP, partZ;
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);
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;
3337 unsigned long long percent = totalelements / 100 *
debug_npercent;
3338 printf(
"total elements = %llu\n", totalelements);
3349 inpart.SetXYZ(1, 0, 0);
3350 for (ir = 0; ir < nrh; ir++)
3352 partR = (ir + 0.5) * deltar + rih;
3355 inpart.SetPerp(partR + deltar);
3357 else if (ir == nrh - 1)
3359 inpart.SetPerp(partR - deltar);
3363 inpart.SetPerp(partR);
3365 for (ip = 0; ip < nph; ip++)
3367 partP = (ip + 0.5) * deltap + pih;
3368 inpart.SetPhi(partP);
3370 for (iz = 0; iz < nzh; iz++)
3372 partZ = (iz) *deltaz + zih;
3375 inpart.SetZ(partZ + deltaz);
3377 else if (iz == nzh - 1)
3379 inpart.SetZ(partZ - deltaz);
3385 partZ += 0.5 * deltaz;
3391 if (
hasTwin && inpart.Z() < 0)
3397 distort =
GetTotalDistortion(inpart.Z() + deltaz, inpart, nSteps,
true, &validToStep);
3399 distort.RotateZ(-inpart.Phi());
3400 diffdistP = distort.Y();
3401 diffdistR = distort.X();
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);
3409 if (
hasTwin && makeUnifiedMap && inpart.Z() < 0)
3417 distortX = distort.X();
3418 distortY = distort.Y();
3419 distort.RotateZ(-inpart.Phi());
3420 distortP = distort.Y();
3421 distortR = distort.X();
3422 distortZ = distort.Z();
3428 hIntDistortionR->Fill(partP, partR, partZ, distortR);
3429 hIntDistortionP->Fill(partP, partR, partZ, distortP);
3430 hIntDistortionZ->Fill(partP, partR, partZ, distortZ);
3434 hIntDistortionX->Fill(partP, partR, partZ, distortX);
3435 hIntDistortionY->Fill(partP, partR, partZ, distortY);
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);
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);
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);
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);
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);
3490 if (!(el % percent))
3493 printf(
"distortion at (ir=%d,ip=%d,iz=%d) is (%E,%E,%E)\n",
3494 ir, ip, iz, distortR, distortP, distortZ);
3501 TCanvas *canvas =
new TCanvas(
"cdistort",
"distortion integrals", 1200, 800);
3504 TPad *
c =
new TPad(
"cplots",
"distortion integral plots", 0, 0.2, 1, 1);
3506 TPad *textpad =
new TPad(
"ctext",
"distortion integral plots", 0, 0.0, 1, 0.2);
3508 gStyle->SetOptStat();
3509 for (
int i = 0; i < 3; i++)
3512 for (
int ax = 0;
ax < 3;
ax++)
3515 c->cd(i * 4 +
ax + 1);
3516 gPad->SetRightMargin(0.15);
3517 hIntDist[
ax][i]->SetStats(0);
3518 hIntDist[
ax][i]->Draw(
"colz");
3521 hRDist[0][i]->SetStats(0);
3522 hRDist[0][i]->SetFillColor(kRed);
3523 hRDist[0][i]->Draw(
"hist");
3526 hRDist[1][i]->SetStats(0);
3527 hRDist[1][i]->SetLineColor(kBlue);
3528 hRDist[1][i]->Draw(
"hist,same");
3547 float texshift = 0.12;
3548 TLatex *tex =
new TLatex(0.0, texpos,
"Fill Me In");
3549 tex->SetTextSize(texshift * 0.8);
3555 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3572 canvas->SaveAs(summaryFilename.Data());
3580 for (
int i = 0; i < 3; i++)
3583 for (
int ax = 0;
ax < 3;
ax++)
3586 c->cd(i * 4 +
ax + 1);
3587 gPad->SetRightMargin(0.15);
3588 hDiffDist[
ax][i]->SetStats(0);
3589 hDiffDist[
ax][i]->Draw(
"colz");
3592 hRDiffDist[0][i]->SetStats(0);
3593 hRDiffDist[0][i]->SetFillColor(kRed);
3594 hRDiffDist[0][i]->Draw(
"hist");
3597 hRDiffDist[1][i]->SetStats(0);
3598 hRDiffDist[1][i]->SetLineColor(kBlue);
3599 hRDiffDist[1][i]->Draw(
"hist same");
3605 tex->SetTextSize(texshift * 0.8);
3611 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3617 tex->DrawLatex(0.05, texpos,
"Differential Plots");
3630 canvas->SaveAs(diffSummaryFilename.Data());
3636 hDistortionR->Write();
3637 hDistortionP->Write();
3638 hDistortionZ->Write();
3639 hIntDistortionR->Write();
3640 hIntDistortionP->Write();
3641 hIntDistortionZ->Write();
3644 hIntDistortionX->Write();
3645 hIntDistortionY->Write();
3672 printf(
"wrote map and summary to %s.\n", filebase);
3673 printf(
"map:%s.\n", distortionFilename.Data());
3674 printf(
"summary:%s.\n", summaryFilename.Data());
3681 int defaultsteps = 100;
3684 return swimToInSteps(zdest, start, defaultsteps, interpolate, &goodtostep);
3697 printf(
"AnnularFieldSim::swimTo asked to swim particle from (%f,%f,%f) which is outside the ROI:\n", start.X(), start.Y(), start.Z());
3699 printf(
"Returning original position.\n");
3705 double zdist = zdest - start.Z();
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);
3723 else if (interpolate)
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());
3743 double EfieldZ = fieldInt.Z() / zdist;
3744 double BfieldZ = fieldIntB.Z() / zdist;
3756 double c0 = 1 / (1 + T2om2);
3757 double c1 = T1om / (1 + T1om * T1om);
3758 double c2 = T2om2 / (1 + T2om2);
3760 TVector3 EintOverEz = 1 / EfieldZ * fieldInt;
3761 TVector3 BintOverBz = 1 / BfieldZ * fieldIntB;
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();
3770 double vprime = (5000 *
cm /
s) / (100 *
V /
cm);
3771 double vdoubleprime = 0;
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());
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);
3788 if (
abs(deltaZ / zdist) > 0.25)
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);
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);
3815 if (!(
abs(deltaX) < 1E3))
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);
3829 TVector3 shift(deltaX, deltaY, deltaZ);
3832 shift.RotateZ(-start.Phi());
3836 shift.RotateZ(start.Phi());
3910 return q->
Get(r, p, z);