20 #include <boost/format.hpp>
28 #define ALMOST_ZERO 0.00001
31 int r,
int roi_r0,
int roi_r1,
int ,
int ,
32 int phi,
int roi_phi0,
int roi_phi1,
int ,
int ,
33 int z,
int roi_z0,
int roi_z1,
int ,
int ,
37 if (roi_r0 >= r || roi_r1 > r || roi_r0 >= roi_r1)
41 if (roi_phi0 >= phi || roi_phi1 > phi || roi_phi0 >= roi_phi1)
43 printf(
"phi roi is out of range or spans the wrap-around. Please spare me that math.\n");
46 if (roi_z0 >= z || roi_z1 > z || roi_z0 >= roi_z1)
54 printf(
"AnnularFieldSim::AnnularFieldSim with (%dx%dx%d) grid\n", r, phi, z);
55 printf(
"units m=%1.2E, cm=%1.2E, mm=%1.2E, um=%1.2E,\n",
m,
cm,
mm,
um);
56 printf(
"units s=%1.2E, us=%1.2E, ns=%1.2E,\n",
s,
us,
ns);
57 printf(
"units C=%1.2E, nC=%1.2E, fC=%1.2E, \n",
C,
nC,
fC);
78 rmin = in_innerRadius *
cm;
79 rmax = in_outerRadius *
cm;
105 printf(
"AnnularFieldSim::AnnularFieldSim set variables nr=%d, nphi=%d, nz=%d\n",
nr,
nphi,
nz);
111 step.SetXYZ(1, 0, 0);
130 printf(
"AnnularFieldSim::AnnularFieldSim set roi variables rmin=%d phimin=%d zmin=%d rmax=%d phimax=%d zmax=%d\n",
169 printf(
"AnnularFieldSim::AnnularFieldSim building Epartial (full3D) with nr_roi=%d nphi_roi=%d nz_roi=%d =~%2.2fM TVector3 objects\n", nr_roi, nphi_roi,
nz_roi,
193 printf(
"lookupCase==HybridRes\n");
203 printf(
"lookupCase==PhiSlice\n");
226 printf(
"lookupCase==Analytic (or NoLookup)\n");
248 printf(
"Ran into wrong lookupCase logic in constructor.\n");
255 int r,
int roi_r0,
int roi_r1,
int in_rLowSpacing,
int in_rHighSize,
256 int phi,
int roi_phi0,
int roi_phi1,
int in_phiLowSpacing,
int in_phiHighSize,
257 int z,
int roi_z0,
int roi_z1,
int in_zLowSpacing,
int in_zHighSize,
260 r, roi_r0, roi_r1, in_rLowSpacing, in_rHighSize,
261 phi, roi_phi0, roi_phi1, in_phiLowSpacing, in_phiHighSize,
262 z, roi_z0, roi_z1, in_zLowSpacing, in_zHighSize,
265 printf(
"AnnularFieldSim::OldConstructor building AnnularFieldSim with ChargeCase::FromFile\n");
269 int r,
int roi_r0,
int roi_r1,
270 int phi,
int roi_phi0,
int roi_phi1,
271 int z,
int roi_z0,
int roi_z1,
274 r, roi_r0, roi_r1, 1, 3,
275 phi, roi_phi0, roi_phi1, 1, 3,
276 z, roi_z0, roi_z1, 1, 3,
281 printf(
"AnnularFieldSim::OldConstructor building AnnularFieldSim with local_size=1 in all dimensions, lowres_spacing=1 in all dimensions\n");
292 printf(
"AnnularFieldSim::OldConstructor building AnnularFieldSim with roi=full in all dimensions\n");
302 printf(
"something's asking for 'at' with r=%f, which is index=%d\n", at.Perp(), r_position);
310 TVector3
field(0, 0, 0);
315 if (
green ==
nullptr)
317 TVector3 delr = at - from;
326 field.SetMag(
k_perm * 1 / (delr * delr));
334 double delphi = abs(at.DeltaPhi(from));
338 double Er =
green->
Er(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
343 Ephi =
green->
Ephi(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
345 double Ez =
green->
Ez(at.Perp(), atphi, at.Z(), from.Perp(), fromphi, from.Z());
346 field.SetXYZ(-Er, -Ephi, -Ez);
348 field.RotateZ(at.Phi());
367 if (p >= 2 * M_PI || p < 0)
369 printf(
"AnnularFieldSim::FilterPhiPos asked to filter %f, which is more than range=%f out of bounds. Check what called this.\n", phi, 2 * M_PI);
391 if (p >= range || p < 0)
393 printf(
"AnnularFieldSim::FilterPhiIndex asked to filter %d, which is more than range=%d out of bounds. Check what called this.\n", phi, range);
401 float r0f = (pos -
rmin) /
step.Perp();
408 float p0f = (
pos) /
step.Phi();
409 int phitemp = floor(p0f);
416 float z0f = (pos -
zmin) /
step.Z();
425 float r0f = (pos -
rmin) /
step.Perp();
451 float p0f = (
pos) /
step.Phi();
452 int phitemp = floor(p0f);
485 float z0f = (pos -
zmin) /
step.Z();
526 if (!rOkay || !phiOkay)
528 printf(
"AnnularFieldSim::analyticFieldIntegral asked to integrate along (r=%f,phi=%f), index=(%d,%d), which is out of bounds. returning starting position.\n", start.Perp(), start.Phi(),
r,
phi);
532 int dir = (start.Z() > zdest ? -1 : 1);
553 bool startOkay = (startBound ==
InBounds);
556 if (!startOkay || !endOkay)
558 printf(
"AnnularFieldSim::analyticFieldIntegral asked to integrate from z=%f to %f, index=%d to %d), which is out of bounds. returning starting position.\n", startz, endz, zi, zf);
585 if (!rOkay || !phiOkay)
587 printf(
"AnnularFieldSim::fieldIntegral asked to integrate along (r=%f,phi=%f), index=(%d,%d), which is out of bounds. returning starting position.\n", start.Perp(), start.Phi(),
r,
phi);
591 int dir = (start.Z() > zdest ? -1 : 1);
612 bool startOkay = (startBound ==
InBounds);
615 if (!startOkay || !endOkay)
617 printf(
"AnnularFieldSim::fieldIntegral asked to integrate from z=%f to %f, index=%d to %d), which is out of bounds. returning starting position.\n", startz, endz, zi, zf);
621 TVector3 fieldInt(0, 0, 0);
624 for (
int i = zi;
i < zf;
i++)
628 fieldInt += tf *
step.Z();
642 return dir * fieldInt;
650 c.SetPerp((r + 0.5) *
step.Perp() +
rmin);
651 c.SetPhi((phi + 0.5) *
step.Phi());
652 c.SetZ((z + 0.5) *
step.Z() +
zmin);
678 float ravg = (r0 +
r1) / 2.0 + 0.5;
685 printf(
"phi1(%d)<=phi0(%d) even after boosting phi1. check what called this!\n", phi1, phi0);
688 float phiavg = (r0 +
r1) / 2.0 + 0.5;
694 float zavg = (z0 + z1) / 2.0 + 0.5;
697 c.SetPerp((ravg) *
step.Perp() +
rmin);
698 c.SetPhi((phiavg) *
step.Phi());
711 float rout = rin +
step.Perp();
713 float rMid = (4 * sin(
step.Phi() / 2) * (pow(rout, 3) - pow(rin, 3)) / (3 *
step.Phi() * (pow(rout, 2) - pow(rin, 2))));
715 c.SetPhi((phi + 0.5) *
step.Phi());
716 c.SetZ((z + 0.5) *
step.Z());
725 float r0 = (start.Perp() -
rmin) /
step.Perp() - 0.5;
727 float r0d = r0 - r0i;
730 ri[2] = ri[3] = r0i + 1;
732 rw[0] = rw[1] = 1 - r0d;
735 bool skip[] = {
false,
false,
false,
false};
750 float p0 = (start.Phi()) /
step.Phi() - 0.5;
752 float p0d = p0 - p0i;
757 pw[0] = pw[2] = 1 - p0d;
774 int dir = (start.Z() < zdest ? 1 : -1);
798 if (!startOkay || !endOkay)
800 printf(
"AnnularFieldSim::InterpolatedFieldIntegral asked to integrate from z=%f to %f, index=%d to %d), which is out of bounds. returning zero\n", startz, endz, zi, zf);
812 bool isE = (field ==
Efield);
813 bool isB = (field ==
Bfield);
814 printf(
"interpFieldInt: isE=%d, isB=%d, start=(%2.4f,%2.4f,%2.4f) dest=%2.4f\n", isE, isB, start.X(), start.Y(), start.Z(), zdest);
815 printf(
"interpolating fieldInt for %s at r=%f,phi=%f\n", (field ==
Efield) ?
"Efield" :
"Bfield", r0, p0);
816 printf(
"direction=%d, startz=%2.4f, endz=%2.4f, zi=%d, zf=%d\n", dir, startz, endz, zi, zf);
817 for (
int i = 0;
i < 4;
i++)
819 printf(
"skip[%d]=%d\tw[i]=%2.2f, (pw=%2.2f, rw=%2.2f)\n",
i, skip[
i], rw[i] * pw[i], pw[i], rw[i]);
823 TVector3 fieldInt(0, 0, 0), partialInt;
825 for (
int i = 0;
i < 4;
i++)
832 partialInt.SetXYZ(0, 0, 0);
833 for (
int j = zi;
j < zf;
j++)
850 fieldInt += rw[
i] * pw[
i] * partialInt;
853 return dir * fieldInt;
865 double ifc_radius = 83.5 *
cm;
866 double ofc_radius = 254.5 *
cm;
867 double tpc_halfz = 250 *
cm;
870 double totalcharge = 0;
871 double localcharge = 0;
874 for (
int ifr = 0; ifr <
nr; ifr++)
876 for (
int ifphi = 0; ifphi <
nphi; ifphi++)
878 for (
int ifz = 0; ifz <
nz; ifz++)
881 pos = pos * (1.0 /
cm);
885 totalcharge += localcharge;
891 printf(
"AnnularFieldSim::load_analytic_spacecharge: Total charge Q=%E Coulombs\n", totalcharge);
903 for (
int ifr = 0; ifr <
nr; ifr++)
906 for (
int ifphi = 0; ifphi <
nphi; ifphi++)
909 for (
int ifz = 0; ifz <
nz; ifz++)
924 TFile fieldFile(filename.c_str(),
"READ");
926 fieldFile.GetObject(treename.c_str(), fTree);
927 Efieldname =
"E:" + filename +
":" + treename;
931 fTree->SetBranchAddress(
"r", &r);
932 fTree->SetBranchAddress(
"er", &fr);
933 fTree->SetBranchAddress(
"z", &z);
934 fTree->SetBranchAddress(
"ez", &fz);
948 TFile fieldFile(filename.c_str(),
"READ");
950 fieldFile.GetObject(treename.c_str(), fTree);
951 Bfieldname =
"B:" + filename +
":" + treename;
954 fTree->SetBranchAddress(
"r", &r);
955 fTree->SetBranchAddress(
"br", &fr);
956 fTree->SetBranchAddress(
"z", &z);
957 fTree->SetBranchAddress(
"bz", &fz);
973 TFile fieldFile(filename.c_str(),
"READ");
975 fieldFile.GetObject(treename.c_str(), fTree);
980 fTree->SetBranchAddress(
"rho", &r);
981 fTree->SetBranchAddress(
"brho", &fr);
982 fTree->SetBranchAddress(
"z", &z);
983 fTree->SetBranchAddress(
"bz", &fz);
984 fTree->SetBranchAddress(
"phi", &phi);
985 fTree->SetBranchAddress(
"bphi", &fphi);
995 bool phiSymmetry = (phiptr ==
nullptr);
996 int lowres_factor = 10;
998 TH3F *htEntries =
new TH3F(
"htentries",
"num of entries in the field loading",
nphi, 0, M_PI * 2.0,
nr,
rmin,
rmax,
nz,
zmin,
zmax);
1000 TH3F *htEntriesLow =
new TH3F(
"htentrieslow",
"num of lowres entries in the field loading",
nphi / lowres_factor + 1, 0, M_PI * 2.0,
nr / lowres_factor + 1,
rmin,
rmax,
nz / lowres_factor + 1,
zmin,
zmax);
1002 char axis[] =
"rpz";
1003 for (
int i = 0;
i < 3;
i++)
1005 htSum[
i] =
new TH3F(Form(
"htsum%d",
i), Form(
"sum of %c-axis entries in the field loading", *(axis +
i)),
nphi, 0, M_PI * 2.0,
nr,
rmin,
rmax,
nz,
zmin,
zmax);
1006 htSumLow[
i] =
new TH3F(Form(
"htsumlow%d",
i), Form(
"sum of low %c-axis entries in the field loading", *(axis +
i)),
nphi / lowres_factor + 1, 0, M_PI * 2.0,
nr / lowres_factor + 1,
rmin,
rmax,
nz / lowres_factor + 1,
zmin,
zmax);
1009 int nEntries = source->GetEntries();
1010 for (
int i = 0;
i < nEntries;
i++)
1012 source->GetEntry(
i);
1013 float zval = *zptr * zsign;
1019 htEntries->Fill(*phiptr, *rptr, zval);
1020 htSum[0]->Fill(*phiptr, *rptr, zval, *frptr * fieldunit);
1021 htSum[1]->Fill(*phiptr, *rptr, zval, *fphiptr * fieldunit);
1022 htSum[2]->Fill(*phiptr, *rptr, zval, *fzptr * fieldunit * zsign);
1023 htEntriesLow->Fill(*phiptr, *rptr, zval);
1024 htSumLow[0]->Fill(*phiptr, *rptr, zval, *frptr * fieldunit);
1025 htSumLow[1]->Fill(*phiptr, *rptr, zval, *fphiptr * fieldunit);
1026 htSumLow[2]->Fill(*phiptr, *rptr, zval, *fzptr * fieldunit * zsign);
1030 for (
int j = 0;
j <
nphi;
j++)
1032 htEntries->Fill(
j *
step.Phi(), *rptr, zval);
1033 htSum[0]->Fill(
j *
step.Phi(), *rptr, zval, *frptr * fieldunit);
1034 htSum[1]->Fill(
j *
step.Phi(), *rptr, zval, *fphiptr * fieldunit);
1035 htSum[2]->Fill(
j *
step.Phi(), *rptr, zval, *fzptr * fieldunit * zsign);
1036 htEntriesLow->Fill(
j *
step.Phi(), *rptr, zval);
1037 htSumLow[0]->Fill(
j *
step.Phi(), *rptr, zval, *frptr * fieldunit);
1038 htSumLow[1]->Fill(
j *
step.Phi(), *rptr, zval, *fphiptr * fieldunit);
1039 htSumLow[2]->Fill(
j *
step.Phi(), *rptr, zval, *fzptr * fieldunit * zsign);
1045 for (
int i = 0;
i <
nphi;
i++)
1047 for (
int j = 0;
j <
nr;
j++)
1049 for (
int k = 0;
k <
nz;
k++)
1052 int bin = htEntries->FindBin(
FilterPhiPos(cellcenter.Phi()), cellcenter.Perp(), cellcenter.Z());
1053 TVector3 fieldvec(htSum[0]->GetBinContent(bin), htSum[1]->GetBinContent(bin), htSum[2]->GetBinContent(bin));
1054 fieldvec = fieldvec * (1.0 / htEntries->GetBinContent(bin));
1055 if (htEntries->GetBinContent(bin) < 0.99)
1062 (*field)->Set(
j,
i,
k, fieldvec);
1068 printf(
"found %d empty bins when constructing %s. Filling with lower resolution.\n", nemptybins, (*field ==
Bfield) ?
"Bfield" :
"Eexternal");
1069 for (
int i = 0;
i <
nphi;
i++)
1071 for (
int j = 0;
j <
nr;
j++)
1073 for (
int k = 0;
k <
nz;
k++)
1076 int bin = htEntries->FindBin(
FilterPhiPos(cellcenter.Phi()), cellcenter.Perp(), cellcenter.Z());
1077 if (htEntries->GetBinContent(bin) == 0)
1079 int lowbin = htEntriesLow->FindBin(
FilterPhiPos(cellcenter.Phi()), cellcenter.Perp(), cellcenter.Z());
1080 TVector3 fieldvec(htSumLow[0]->GetBinContent(lowbin), htSumLow[1]->GetBinContent(lowbin), htSumLow[2]->GetBinContent(lowbin));
1081 fieldvec = fieldvec * (1.0 / htEntriesLow->GetBinContent(lowbin));
1082 if (htEntriesLow->GetBinContent(lowbin) < 0.99)
1084 printf(
"not enough entries in source to fill fieldmap. None near r=%f, phi=%f, z=%f. Pick lower granularity!\n",
1085 cellcenter.Perp(),
FilterPhiPos(cellcenter.Phi()), cellcenter.Z());
1090 (*field)->Set(
j,
i,
k, fieldvec);
1101 TFile *
f = TFile::Open(filename.c_str());
1103 f->GetObject(histname.c_str(), scmap);
1104 std::cout <<
"Loading spacecharge from '" << filename
1105 <<
"'. Seeking histname '" << histname <<
"'" << std::endl;
1115 TFile *
f = TFile::Open(filename.c_str());
1117 f->GetObject(histname.c_str(), scmap);
1118 std::cout <<
"Resampling spacecharge from '" << filename
1119 <<
"'. Seeking histname '" << histname <<
"'" << std::endl;
1135 float hrmin = hist->GetYaxis()->GetXmin();
1136 float hrmax = hist->GetYaxis()->GetXmax();
1137 float hphimin = hist->GetXaxis()->GetXmin();
1138 float hphimax = hist->GetXaxis()->GetXmax();
1139 float hzmin = hist->GetZaxis()->GetXmin();
1140 float hzmax = hist->GetZaxis()->GetXmax();
1143 int hrn = hist->GetNbinsY();
1144 int hphin = hist->GetNbinsX();
1145 int hzn = hist->GetNbinsZ();
1146 printf(
"From histogram, native r dimensions: %f<r<%f, hrn (Y axis)=%d\n", hrmin, hrmax, hrn);
1147 printf(
"From histogram, native phi dimensions: %f<phi<%f, hrn (X axis)=%d\n", hphimin, hphimax, hphin);
1148 printf(
"From histogram, native z dimensions: %f<z<%f, hrn (Z axis)=%d\n", hzmin, hzmax, hzn);
1151 float hrstep = (hrmax - hrmin) / hrn;
1152 float hphistep = (hphimax - hphimin) / hphin;
1153 float hzstep = (hzmax - hzmin) / hzn;
1154 float halfrstep = hrstep * 0.5;
1155 float halfphistep = hphistep * 0.5;
1156 float halfzstep = hzstep * 0.5;
1161 if (!isChargeDensity)
1163 for (
int i = 0;
i < hphin;
i++)
1165 float phi = hphimin + hphistep *
i;
1166 for (
int j = 0;
j < hrn;
j++)
1168 float r = hrmin + hrstep *
j;
1169 for (
int k = 0;
k < hzn;
k++)
1171 float z = hzmin + hzstep *
k;
1172 int bin = hist->FindBin(phi + halfphistep, r + halfrstep, z + halfzstep);
1173 double vol = hzstep * hphistep * (r + halfrstep) * hrstep;
1174 hist->SetBinContent(bin, hist->GetBinContent(bin) / vol);
1180 TH3F *resampled =
new TH3F(
"resampled",
"resampled charge", new_nphi, hphimin, hphimax, new_nr, hrmin, hrmax, new_nz, hzmin, hzmax);
1181 float new_phistep = (hphimax - hphimin) / new_nphi;
1182 float new_rstep = (hrmax - hrmin) / new_nr;
1183 float new_zstep = (hzmax - hzmin) / new_nz;
1186 for (
int i = 0;
i < new_nphi;
i++)
1188 float phi = hphimin + new_phistep * (
i + 0.5);
1189 float hphi = (phi - hphimin) / hphistep;
1190 bool phisafe = ((hphi - 0.75) > 0) && ((hphi + 0.75) < hphin);
1191 for (
int j = 0;
j < new_nr;
j++)
1193 float r = hrmin + new_rstep * (
j + 0.5);
1194 float hr = (r - hrmin) / hrstep;
1195 bool rsafe = ((hr - 0.5) > 0) && ((hr + 0.5) < hrn);
1196 for (
int k = 0;
k < new_nz;
k++)
1198 float z = hzmin + new_zstep * (
k + 0.5);
1199 float hz = (z - hzmin) / hzstep;
1200 bool zsafe = ((hz - 0.5) > 0) && ((hz + 0.5) < hzn);
1202 if (phisafe && rsafe && zsafe)
1205 resampled->Fill(phi, r, z, hist->Interpolate(phi, r, z));
1209 int bin = hist->FindBin(phi, r, z);
1210 resampled->Fill(phi, r, z, hist->GetBinContent(bin));
1215 printf(
"here in load_and_resample_spacecharge, I am going to call load_spacecharge.\n");
1222 if (abs(zoffset) > 0.001)
1224 printf(
"nonzero zoffset given (%E) but new spacecharge loader can't deal with that. Failing.\n", zoffset);
1227 if (isChargeDensity)
1229 printf(
"Input dataset is flagged as recording density not total charge, but new loader can't deal with that Failing..\n");
1234 sprintf(
chargestring,
"SC loaded externally: %s.", inputchargestring);
1242 sprintf(
chargestring,
"SC loaded externally: %s.", inputchargestring);
1249 TH3F *hsc =
new TH3F(
"hInternalSpaceCharge",
"Internal Charge Histogram;phi(rad);r(cm);z(cm)",
nphi, 0,
phispan,
nr,
rmin,
rmax,
nz,
zmin,
zmax);
1250 for (
int i = 0;
i <
nphi;
i++)
1252 float phi = 0 +
step.Phi() * (
i + 0.5);
1253 for (
int j = 0;
j <
nr;
j++)
1256 for (
int k = 0;
k <
nz;
k++)
1264 hsc->SaveAs(filename.c_str());
1333 unsigned long long totalelements =
nr_roi;
1336 unsigned long long percent = totalelements / 100 *
debug_npercent;
1337 printf(
"total elements = %llu\n", totalelements *
nr *
nphi *
nz);
1349 if (!(el % percent))
1352 printf(
"sum_field_at (ir=%d,iphi=%d,iz=%d) gives (%E,%E,%E)\n",
1353 ir, iphi, iz, localF.X(), localF.Y(), localF.Z());
1375 printf(
"lookupCase==Full3D\n");
1381 printf(
"lookupCase==HybridRes\n");
1387 printf(
"Populating lookup: lookupCase==PhiSlice\n");
1392 printf(
"Populating lookup: lookupCase==Analytic ===> skipping!\n");
1396 printf(
"Populating lookup: lookupCase==NoLookup ===> skipping!\n");
1410 printf(
"populating full lookup table for (%dx%dx%d)x(%dx%dx%d) grid\n",
1415 totalelements *=
nr;
1416 totalelements *=
nphi;
1417 totalelements *=
nz;
1418 unsigned long long percent = totalelements / 100 *
debug_npercent;
1419 printf(
"total elements = %llu\n", totalelements);
1420 TVector3
at(1, 0, 0);
1421 TVector3 from(1, 0, 0);
1422 TVector3
zero(0, 0, 0);
1432 for (
int ior = 0; ior <
nr; ior++)
1434 for (
int iophi = 0; iophi <
nphi; iophi++)
1436 for (
int ioz = 0; ioz <
nz; ioz++)
1439 if (!(el % percent))
1447 if (ifr == ior && ifphi == iophi && ifz == ioz)
1466 TVector3
at(1, 0, 0);
1467 TVector3 from(1, 0, 0);
1468 TVector3
zero(0, 0, 0);
1471 int r_highres_dist = (
nr_high - 1) / 2;
1472 int phi_highres_dist = (
nphi_high - 1) / 2;
1473 int z_highres_dist = (
nz_high - 1) / 2;
1476 static int nfbinsin[3][3][3];
1477 for (
auto &
i : nfbinsin)
1479 for (
int j = 0;
j < 3;
j++)
1481 for (
int k = 0;
k < 3;
k++)
1491 TVector3 currentf, newf;
1494 int r_parentlow = floor((ifr - r_highres_dist) / (
r_spacing * 1.0));
1495 int r_parenthigh = floor((ifr + r_highres_dist) / (
r_spacing * 1.0)) + 1;
1496 int r_startpoint = r_parentlow *
r_spacing;
1497 int r_endpoint = r_parenthigh *
r_spacing;
1502 bool phi_parentlow_wrapped = (ifphi - phi_highres_dist < 0);
1504 if (phi_parentlow_wrapped)
1506 phi_startpoint -=
nphi;
1509 int phi_parenthigh = floor(
FilterPhiIndex(ifphi + phi_highres_dist) / (phi_spacing * 1.0)) + 1;
1510 bool phi_parenthigh_wrapped = (ifphi + phi_highres_dist >=
nphi);
1512 if (phi_parenthigh_wrapped)
1514 phi_endpoint +=
nphi;
1521 int z_parentlow = floor((ifz - z_highres_dist) / (
z_spacing * 1.0));
1522 int z_parenthigh = floor((ifz + z_highres_dist) / (
z_spacing * 1.0)) + 1;
1523 int z_startpoint = z_parentlow *
z_spacing;
1524 int z_endpoint = z_parenthigh *
z_spacing;
1535 for (
int ir = r_startpoint; ir < r_endpoint; ir++)
1548 int rbin = (ir - ifr) + r_highres_dist;
1561 for (
int iphi = phi_startpoint; iphi < phi_endpoint; iphi++)
1565 int phibin = (iphi - ifphi) + phi_highres_dist;
1577 for (
int iz = z_startpoint; iz < z_endpoint; iz++)
1587 int zbin = (iz - ifz) + z_highres_dist;
1602 nfbinsin[rcell][phicell][zcell]++;
1603 int nf = nfbinsin[rcell][phicell][zcell];
1608 if (zcell != 1 || rcell != 1 || phicell != 1)
1614 printf(
"%d: Getting with phi=%d\n", __LINE__, iphi_rel);
1619 newf = (currentf * (nf - 1) +
calc_unit_field(at, from)) * (1 / (nf * 1.0));
1626 if (ifr == rbin && ifphi == phibin && ifz == zbin)
1647 TVector3
at(1, 0, 0);
1648 TVector3 from(1, 0, 0);
1649 TVector3
zero(0, 0, 0);
1650 int fphi_low, fphi_high, fz_low, fz_high;
1651 int r_low, r_high, phi_low, phi_high, z_low, z_high;
1657 int fr_high = fr_low + r_spacing - 1;
1665 fphi_high = fphi_low + phi_spacing - 1;
1666 if (fphi_high >=
nphi)
1668 fphi_high =
nphi - 1;
1673 fz_high = fz_low + z_spacing - 1;
1682 for (
int ior = 0; ior <
nr_low; ior++)
1685 r_high = r_low + r_spacing - 1;
1692 for (
int iophi = 0; iophi <
nphi_low; iophi++)
1695 phi_high = phi_low + phi_spacing - 1;
1696 if (phi_high >=
nphi)
1698 phi_high =
nphi - 1;
1701 for (
int ioz = 0; ioz <
nz_low; ioz++)
1704 z_high = z_low + z_spacing - 1;
1712 if (ifr == ior && ifphi == iophi && ifz == ioz)
1739 unsigned long long totalelements =
nr;
1740 totalelements *=
nphi;
1741 totalelements *=
nz;
1744 unsigned long long percent = totalelements / 100 *
debug_npercent;
1745 printf(
"total elements = %llu\n", totalelements);
1746 TVector3
at(1, 0, 0);
1747 TVector3 from(1, 0, 0);
1748 TVector3
zero(0, 0, 0);
1756 for (
int ior = 0; ior <
nr; ior++)
1758 for (
int iophi = 0; iophi <
nphi; iophi++)
1760 for (
int ioz = 0; ioz <
nz; ioz++)
1766 if (ifr == ior && 0 == iophi && ifz == ioz)
1768 if (!(el % percent))
1771 printf(
"self-to-self is zero (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) gives (%E,%E,%E)\n",
1772 ior, iophi, ioz, ifr, ifz, zero.X(), zero.Y(), zero.Z());
1781 if (!(el % percent))
1784 printf(
"calc_unit_field (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) gives (%E,%E,%E)\n",
1785 ior, iophi, ioz, ifr, ifz, unitf.X(), unitf.Y(), unitf.Z());
1802 unsigned long long totalelements =
nr;
1803 totalelements *=
nphi;
1804 totalelements *=
nz;
1807 unsigned long long percent = totalelements / 100 *
debug_npercent;
1808 printf(
"total elements = %llu\n", totalelements);
1810 TFile *
input = TFile::Open(sourcefile,
"READ");
1812 input->GetObject(
"info", tInfo);
1815 float file_rmin, file_rmax, file_zmin, file_zmax;
1816 int file_rmin_roi, file_rmax_roi, file_zmin_roi, file_zmax_roi;
1817 int file_nr, file_np, file_nz;
1818 tInfo->SetBranchAddress(
"rmin", &file_rmin);
1819 tInfo->SetBranchAddress(
"rmax", &file_rmax);
1820 tInfo->SetBranchAddress(
"zmin", &file_zmin);
1821 tInfo->SetBranchAddress(
"zmax", &file_zmax);
1822 tInfo->SetBranchAddress(
"rmin_roi_index", &file_rmin_roi);
1823 tInfo->SetBranchAddress(
"rmax_roi_index", &file_rmax_roi);
1824 tInfo->SetBranchAddress(
"zmin_roi_index", &file_zmin_roi);
1825 tInfo->SetBranchAddress(
"zmax_roi_index", &file_zmax_roi);
1826 tInfo->SetBranchAddress(
"nr", &file_nr);
1827 tInfo->SetBranchAddress(
"nphi", &file_np);
1828 tInfo->SetBranchAddress(
"nz", &file_nz);
1831 printf(
"param\tobj\tfile\n");
1832 printf(
"nr\t%d\t%d\n",
nr, file_nr);
1834 printf(
"nz\t%d\t%d\n",
nz, file_nz);
1835 printf(
"rmin\t%2.2f\t%2.2f\n",
rmin, file_rmin);
1836 printf(
"rmax\t%2.2f\t%2.2f\n",
rmax, file_rmax);
1837 printf(
"zmin\t%2.2f\t%2.2f\n",
zmin, file_zmin);
1838 printf(
"zmax\t%2.2f\t%2.2f\n",
zmax, file_zmax);
1844 if (file_rmin !=
rmin || file_rmax !=
rmax ||
1845 file_zmin !=
zmin || file_zmax !=
zmax ||
1848 file_nr !=
nr || file_np !=
nphi || file_nz !=
nz)
1850 printf(
"file parameters do not match fieldsim parameters:\n");
1856 input->GetObject(
"phislice", tLookup);
1858 int ior, ifr, iophi, ioz, ifz;
1859 TVector3 *unitf =
nullptr;
1860 tLookup->SetBranchAddress(
"ir_source", &ior);
1861 tLookup->SetBranchAddress(
"ir_target", &ifr);
1862 tLookup->SetBranchAddress(
"ip_source", &iophi);
1864 tLookup->SetBranchAddress(
"iz_source", &ioz);
1865 tLookup->SetBranchAddress(
"iz_target", &ifz);
1866 tLookup->SetBranchAddress(
"Evec", &unitf);
1869 printf(
"%s has %lld entries\n", sourcefile, tLookup->GetEntries());
1870 for (
int i = 0;
i < (int) totalelements;
i++)
1873 tLookup->GetEntry(
i);
1877 if (!(el % percent))
1880 printf(
"field from (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) is (%E,%E,%E)\n",
1881 ior, iophi, ioz, ifr, ifz, unitf->X(), unitf->Y(), unitf->Z());
1892 unsigned long long totalelements =
nr;
1893 totalelements *=
nphi;
1894 totalelements *=
nz;
1897 unsigned long long percent = totalelements / 100 *
debug_npercent;
1898 printf(
"total elements = %llu\n", totalelements);
1900 TFile *
output = TFile::Open(destfile,
"RECREATE");
1903 TTree *tInfo =
new TTree(
"info",
"Information about Lookup Table");
1904 tInfo->Branch(
"rmin", &
rmin);
1905 tInfo->Branch(
"rmax", &
rmax);
1906 tInfo->Branch(
"zmin", &
zmin);
1907 tInfo->Branch(
"zmax", &
zmax);
1908 tInfo->Branch(
"rmin_roi_index", &
rmin_roi);
1909 tInfo->Branch(
"rmax_roi_index", &
rmax_roi);
1910 tInfo->Branch(
"zmin_roi_index", &
zmin_roi);
1911 tInfo->Branch(
"zmax_roi_index", &
zmax_roi);
1912 tInfo->Branch(
"nr", &
nr);
1913 tInfo->Branch(
"nphi", &
nphi);
1914 tInfo->Branch(
"nz", &
nz);
1916 printf(
"info tree built.\n");
1918 TTree *tLookup =
new TTree(
"phislice",
"Phislice Lookup Table");
1919 int ior, ifr, iophi, ioz, ifz;
1921 tLookup->Branch(
"ir_source", &ior);
1922 tLookup->Branch(
"ir_target", &ifr);
1923 tLookup->Branch(
"ip_source", &iophi);
1925 tLookup->Branch(
"iz_source", &ioz);
1926 tLookup->Branch(
"iz_target", &ifz);
1927 tLookup->Branch(
"Evec", &unitf);
1928 printf(
"lookup tree built.\n");
1935 for (ior = 0; ior <
nr; ior++)
1937 for (iophi = 0; iophi <
nphi; iophi++)
1939 for (ioz = 0; ioz <
nz; ioz++)
1945 if (!(el % percent))
1948 printf(
"field from (ir=%d,iphi=%d,iz=%d) to (or=%d,ophi=0,oz=%d) is (%E,%E,%E)\n",
1949 ior, iophi, ioz, ifr, ifz, unitf.X(), unitf.Y(), unitf.Z());
1970 printf(
"AnnularFieldSim::setFlatFields(B=%f Tesla,E=%f V/cm)\n", B, E);
1973 sprintf(fieldstr,
"%f", E);
1975 sprintf(fieldstr,
"%f", B);
2006 TVector3
sum(0, 0, 0);
2031 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);
2055 for (
int iphi = 0; iphi <
nphi; iphi++)
2059 phidist = fmin(abs(iphi - phi), abs(abs(iphi - phi) - nphi));
2066 for (
int iz = 0; iz <
nz; iz++)
2070 zdist = abs(iz - z);
2078 if (r == ir && phi == iphi && z == iz)
2111 int r_highres_dist = (
nr_high - 1) / 2;
2112 int phi_highres_dist = (
nphi_high - 1) / 2;
2113 int z_highres_dist = (
nz_high - 1) / 2;
2115 int r_parentlow = floor((r - r_highres_dist) / (
r_spacing * 1.0));
2116 int r_parenthigh = floor((r + r_highres_dist) / (
r_spacing * 1.0)) + 1;
2117 int phi_parentlow = floor((phi - phi_highres_dist) / (
phi_spacing * 1.0));
2118 int phi_parenthigh = floor((phi + phi_highres_dist) / (
phi_spacing * 1.0)) + 1;
2119 int z_parentlow = floor((z - z_highres_dist) / (
z_spacing * 1.0));
2120 int z_parenthigh = floor((z + z_highres_dist) / (
z_spacing * 1.0)) + 1;
2121 printf(
"AnnularFieldSim::sum_local_field_at parents: rlow=%d,philow=%d,zlow=%d,rhigh=%d,phihigh=%d,zhigh=%d\n", r_parentlow, phi_parentlow, z_parentlow, r_parenthigh, phi_parenthigh, z_parenthigh);
2140 int rbin = (ir -
r) + r_highres_dist;
2153 int phibin = (iphi -
phi) + phi_highres_dist;
2174 int zbin = (iz -
z) + z_highres_dist;
2194 TVector3
sum(0, 0, 0);
2202 for (
int ir = 0; ir <
nr_high; ir++)
2204 for (
int iphi = 0; iphi <
nphi_high; iphi++)
2206 for (
int iz = 0; iz <
nz_high; iz++)
2228 bool skip[] = {
false,
false,
false,
false,
false,
false,
false,
false};
2231 int r0i = floor(r0);
2232 float r0d = r0 - r0i;
2244 for (
int i = 0;
i < 8;
i++)
2246 if ((
i / 4) % 2 == 0)
2255 for (
int i = 0;
i < 8;
i++)
2257 if ((
i / 4) % 2 == 1)
2270 int p0i = floor(p0);
2271 float p0d = p0 - p0i;
2290 for (
int i = 0;
i < 8;
i++)
2292 if ((
i / 2) % 2 == 0)
2301 for (
int i = 0;
i < 8;
i++)
2303 if ((
i / 2) % 2 == 1)
2314 int z0i = floor(z0);
2315 float z0d = z0 - z0i;
2326 for (
int i = 0;
i < 8;
i++)
2337 for (
int i = 0;
i < 8;
i++)
2347 TVector3
sum(0, 0, 0);
2352 bool overlapsPhi, overlapsZ;
2354 int r_highres_dist = (
nr_high - 1) / 2;
2355 int phi_highres_dist = (
nphi_high - 1) / 2;
2356 int z_highres_dist = (
nz_high - 1) / 2;
2358 for (
int ir = 0; ir <
nr_low; ir++)
2361 lBinEdge[1] = (ir + 1) * r_spacing - 1;
2362 hRegionEdge[0] = r - r_highres_dist;
2363 hRegionEdge[1] = r + r_highres_dist;
2364 bool overlapsR = (lBinEdge[0] <= hRegionEdge[1] && hRegionEdge[0] <= lBinEdge[1]);
2365 for (
int iphi = 0; iphi <
nphi_low; iphi++)
2368 lBinEdge[1] = (iphi + 1) * phi_spacing - 1;
2369 hRegionEdge[0] = phi - phi_highres_dist;
2370 hRegionEdge[1] = phi + phi_highres_dist;
2371 overlapsPhi = (lBinEdge[0] <= hRegionEdge[1] && hRegionEdge[0] <= lBinEdge[1]);
2372 for (
int iz = 0; iz <
nz_low; iz++)
2375 lBinEdge[1] = (iz + 1) * z_spacing - 1;
2376 hRegionEdge[0] = z - z_highres_dist;
2377 hRegionEdge[1] = z + z_highres_dist;
2378 overlapsZ = (lBinEdge[0] <= hRegionEdge[1] && hRegionEdge[0] <= lBinEdge[1]);
2382 if (overlapsR && overlapsPhi && overlapsZ)
2391 for (
int i = 0;
i < 8;
i++)
2399 printf(
"considering an l-bins effect on itself, r=%d,phi=%d,z=%d (matches i=%d, not skipped), means we're not interpolating fairly\n", ir, iphi, iz, i);
2404 if (pi[(i / 2) % 2] < 0)
2406 printf(
"%d: Getting with phi=%d\n", __LINE__, pi[(i / 2) % 2]);
2408 sum += (
Epartial_lowres->
Get(ri[(i / 4) % 2], pi[(i / 2) % 2], zi[(i) % 2], ir, iphi, iz) *
q_lowres->
Get(ir, iphi, iz)) * zw[(i) % 2] * pw[(i / 2) % 2] * rw[(i / 4) % 2];
2425 float rotphi = pos.Phi() - slicepos.Phi();
2433 TVector3
sum(0, 0, 0);
2434 TVector3 unrotatedField(0, 0, 0);
2435 TVector3 unitField(0, 0, 0);
2437 for (
int ir = 0; ir <
nr; ir++)
2439 for (
int iphi = 0; iphi <
nphi; iphi++)
2441 for (
int iz = 0; iz <
nz; iz++)
2444 if (r == ir && phi == iphi && z == iz)
2450 unitField.RotateZ(rotphi);
2472 double zdist = (zdest - start.Z()) *
cm;
2473 double zstep = zdist /
steps;
2476 TVector3 ret =
start;
2477 TVector3 accumulated_distortion(0, 0, 0);
2478 TVector3 accumulated_drift(0, 0, 0);
2479 TVector3 drift_step(0, 0, zstep);
2494 "AnnularFieldSim::swimToInAnalyticSteps at step %d,"
2495 "asked to swim particle from (%f,%f,%f)(cm) (rphiz)=(%fcm,%frad,%fcm)which is outside the ROI.\n",
2496 i, ret.X() /
cm, ret.Y() /
cm, ret.Z() /
cm, ret.Perp() /
cm, ret.Phi(), ret.Z() /
cm);
2497 printf(
" -- %f <= r < %f \t%f <= phi < %f \t%f <= z < %f \n",
2501 printf(
"Returning last valid position.\n");
2502 if (!(goodToStep ==
nullptr))
2504 *goodToStep =
i - 1;
2506 return ret * (1.0 /
cm);
2512 accumulated_drift += drift_step;
2515 ret = start + accumulated_distortion + accumulated_drift;
2518 return ret * (1.0 /
cm);
2523 int *success =
nullptr;
2524 TVector3 straightline(start.X(), start.Y(), zdest);
2526 return straightline + distortion;
2553 printf(
"AnnularFieldSim::GetTotalDistortion starting at (%f,%f,%f)=(r%f,p%f,z%f) asked to drift to z=%f, which is outside the ROI. hasTwin= %d. Returning zero_vector.\n", start.X(), start.Y(), start.Z(), start.Perp(),
FilterPhiPos(start.Phi()), start.Z(), zdest, (int)
hasTwin);
2567 printf(
"AnnularFieldSim::GetTotalDistortion starting at (%f,%f,%f)=(r%f,p%f,z%f) asked to drift from z=%f, which is outside the ROI. Returning zero_vector.\n", start.X(), start.Y(), start.Z(), start.Perp(),
FilterPhiPos(start.Phi()), start.Z(), start.Z());
2581 double zstep = (
zmax -
zmin) / steps;
2586 if ((zdest - start.Z()) < 0)
2590 int integernumbersteps = (zdest - start.Z()) / zstep;
2602 TVector3 accumulated_distortion(0, 0, 0);
2603 TVector3 accumulated_drift(0, 0, 0);
2604 TVector3 drift_step(0, 0, zstep);
2612 accumulated_distortion +=
GetStepDistortion(zdest - (integernumbersteps * zstep), position,
true,
false);
2616 position.SetX(start.X() + accumulated_distortion.X());
2617 position.SetY(start.Y() + accumulated_distortion.Y());
2618 position.SetZ(zdest - (integernumbersteps * zstep));
2620 float StartR, FinalR;
2622 for (
int i = 0;
i < integernumbersteps;
i++)
2629 if (!(goodToStep ==
nullptr))
2631 *goodToStep =
i - 1, *success = 0;
2634 return (accumulated_distortion);
2637 StartR = position.Perp();
2640 position.SetX(start.X() + accumulated_distortion.X());
2641 position.SetY(start.Y() + accumulated_distortion.Y());
2642 FinalR = position.Perp();
2647 DeltaR = FinalR - StartR;
2652 position.SetZ(position.Z() + zstep);
2656 *goodToStep = integernumbersteps - 1;
2658 return (accumulated_distortion);
2660 *goodToStep = integernumbersteps;
2662 return accumulated_distortion;
2667 bool mapEfield =
true;
2672 which = mapEfield ?
'E' :
'B';
2676 sprintf(units,
"V/cm");
2680 sprintf(units,
"T");
2683 printf(
"plotting field slices for %c field, slicing at (%1.2F,%1.2f,%1.2f)...\n", which, pos.Perp(),
FilterPhiPos(pos.Phi()), pos.Z());
2684 std::cout <<
"file=" << filebase << std::endl;
2686 TString plotfilename = TString::Format(
"%s.%cfield_slices.pdf", filebase, which);
2693 TH2F *hEfield[3][3];
2695 TH1F *hEfieldComp[3][3];
2696 char axis[] =
"rpzrpz";
2697 float axisval[] = {(float) pos.Perp(), (float)
FilterPhiPos(pos.Phi()), (
float) pos.Z(), (float) pos.Perp(), (float)
FilterPhiPos(pos.Phi()), (
float) pos.Z()};
2699 float axtop[] = {(float) outer.Perp(), 2 * M_PI, (float) outer.Z(), (float) outer.Perp(), 2 * M_PI, (float) outer.Z()};
2700 float axbot[] = {(float) inner.Perp(), 0, (float) inner.Z(), (float) inner.Perp(), 0, (float) inner.Z()};
2708 printf(
"rpz bounds are %f<r%f\t %f<phi%f\t %f<z%f\n", axbot[0], axtop[0], axbot[1], axtop[1], axbot[2], axtop[2]);
2710 for (
int i = 0;
i < 6;
i++)
2712 axstep[
i] = (axtop[
i] - axbot[
i]) / (1.0 * axn[
i]);
2718 for (
int ax = 0;
ax < 3;
ax++)
2721 hCharge[
ax] =
new TH2F(Form(
"hCharge%c", axis[
ax]),
2722 Form(
"Spacecharge Distribution in the %c%c plane at %c=%2.3f (C/cm^3);%c;%c",
2723 axis[ax + 1], axis[ax + 2], axis[ax], axisval[ax], axis[ax + 1], axis[ax + 2]),
2724 axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2725 axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2726 for (
int i = 0;
i < 3;
i++)
2729 hEfield[
ax][
i] =
new TH2F(Form(
"h%cfield%c_%c%c", which, axis[
i], axis[ax + 1], axis[ax + 2]),
2730 Form(
"%c component of %c Field in the %c%c plane at %c=%2.3f (%s);%c;%c",
2731 axis[i], which, axis[ax + 1], axis[ax + 2], axis[ax], axisval[ax], units, axis[ax + 1], axis[ax + 2]),
2732 axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2733 axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2734 hEfieldComp[
ax][
i] =
new TH1F(Form(
"h%cfieldComp%c_%c%c", which, axis[i], axis[ax + 1], axis[ax + 2]),
2735 Form(
"Log Magnitude of %c component of %c Field in the %c%c plane at %c=%2.3f;log10(mag)",
2736 axis[i], which, axis[ax + 1], axis[ax + 2], axis[ax], axisval[ax]),
2742 for (
int ax = 0;
ax < 3;
ax++)
2744 rpz_coord[
ax] = axisval[
ax] + axstep[
ax] / 2;
2745 for (
int i = 0;
i < axn[
ax + 1];
i++)
2747 rpz_coord[(
ax + 1) % 3] = axbot[
ax + 1] + (
i + 0.5) * axstep[
ax + 1];
2748 for (
int j = 0;
j < axn[
ax + 2];
j++)
2750 rpz_coord[(
ax + 2) % 3] = axbot[
ax + 2] + (
j + 0.5) * axstep[
ax + 2];
2751 lpos.SetXYZ(rpz_coord[0], 0, rpz_coord[2]);
2752 lpos.SetPhi(rpz_coord[1]);
2753 if (
false &&
ax == 0)
2755 printf(
"sampling rpz=(%f,%f,%f)=(%f,%f,%f) after conversion to xyz=(%f,%f,%f)\n",
2756 rpz_coord[0], rpz_coord[1], rpz_coord[2],
2757 lpos.Perp(),
FilterPhiPos(lpos.Phi()), lpos.Z(), lpos.X(), lpos.Y(), lpos.Z());
2769 field.RotateZ(-rpz_coord[1]);
2771 hEfield[
ax][0]->Fill(rpz_coord[(
ax + 1) % 3], rpz_coord[(
ax + 2) % 3], field.X());
2772 hEfield[
ax][1]->Fill(rpz_coord[(
ax + 1) % 3], rpz_coord[(
ax + 2) % 3], field.Y());
2773 hEfield[
ax][2]->Fill(rpz_coord[(
ax + 1) % 3], rpz_coord[(
ax + 2) % 3], field.Z());
2774 hCharge[
ax]->Fill(rpz_coord[(
ax + 1) % 3], rpz_coord[(
ax + 2) % 3],
GetChargeAt(lpos));
2775 hEfieldComp[
ax][0]->Fill((abs(field.X())));
2776 hEfieldComp[
ax][1]->Fill((abs(field.Y())));
2777 hEfieldComp[
ax][2]->Fill((abs(field.Z())));
2782 c =
new TCanvas(
"cfieldslices",
"electric field", 1200, 800);
2784 gStyle->SetOptStat();
2785 for (
int ax = 0;
ax < 3;
ax++)
2787 for (
int i = 0;
i < 3;
i++)
2789 c->cd(
ax * 4 +
i + 1);
2790 gPad->SetRightMargin(0.15);
2791 hEfield[
ax][
i]->SetStats(
false);
2792 hEfield[
ax][
i]->Draw(
"colz");
2798 gPad->SetRightMargin(0.15);
2799 hCharge[
ax]->SetStats(
false);
2800 hCharge[
ax]->Draw(
"colz");
2809 c->SaveAs(plotfilename);
2810 printf(
"after plotting field slices...\n");
2811 std::cout <<
"file=" << filebase << std::endl;
2822 TVector3 steplocal(
step.Perp() / r_subsamples, 0,
step.Z() / z_subsamples);
2823 steplocal.SetPhi(
step.Phi() / p_subsamples);
2824 float deltar = steplocal.Perp();
2825 float deltap = steplocal.Phi();
2826 float deltaz = steplocal.Z();
2828 TVector3 stepzvec(0, 0, deltaz);
2848 int nph =
nphi * p_subsamples + 2;
2849 int nrh =
nr * r_subsamples + 2;
2850 int nzh =
nz * z_subsamples + 2;
2852 float rih = lowerEdge.Perp() - 0.5 *
step.Perp() - steplocal.Perp();
2853 float rfh = upperEdge.Perp() + 0.5 *
step.Perp() + steplocal.Perp();
2854 float pih =
FilterPhiPos(lowerEdge.Phi()) - 0.5 *
step.Phi() - steplocal.Phi();
2855 float pfh =
FilterPhiPos(upperEdge.Phi()) + 0.5 *
step.Phi() + steplocal.Phi();
2856 float zih = lowerEdge.Z() - 0.5 *
step.Z() - steplocal.Z();
2857 float zfh = upperEdge.Z() + 0.5 *
step.Z() + steplocal.Z();
2858 float z_readout = upperEdge.Z() + 0.5 *
step.Z();
2860 printf(
"generating distortion map...\n");
2861 printf(
"file=%s\n", filebase);
2865 TString distortionFilename;
2866 distortionFilename.Form(
"%s.distortion_map.hist.root", filebase);
2867 TString summaryFilename;
2868 summaryFilename.Form(
"%s.distortion_summary.pdf", filebase);
2869 TString diffSummaryFilename;
2870 diffSummaryFilename.Form(
"%s.differential_summary.pdf", filebase);
2872 TFile *outf = TFile::Open(distortionFilename.Data(),
"RECREATE");
2877 TH3F *hDistortionR =
new TH3F(
"hDistortionR",
"Per-z-bin Distortion in the R direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2878 TH3F *hDistortionP =
new TH3F(
"hDistortionP",
"Per-z-bin Distortion in the RPhi direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2879 TH3F *hDistortionZ =
new TH3F(
"hDistortionZ",
"Per-z-bin Distortion in the Z direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2880 TH3F *hIntDistortionR =
new TH3F(
"hIntDistortionR",
"Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2881 TH3F *hIntDistortionP =
new TH3F(
"hIntDistortionP",
"Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2882 TH3F *hIntDistortionZ =
new TH3F(
"hIntDistortionZ",
"Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2884 TH3F *hIntDistortionX =
new TH3F(
"hIntDistortionX",
"Integrated X Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2885 TH3F *hIntDistortionY =
new TH3F(
"hIntDistortionY",
"Integrated Y Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
2887 const int nMapComponents = 6;
2888 TH3F *hSeparatedMapComponent[2][6];
2889 TH3F *hValidToStepComponent[2];
2890 TH3F *hReachesReadout[2];
2898 TString sepAxis[] = {
"X",
"Y",
"Z",
"R",
"P",
"RPhi"};
2899 float zlower, zupper;
2900 for (
int i = 0;
i < nSides;
i++)
2904 zlower = fmin(zih, zfh);
2905 zupper = fmax(zih, zfh);
2909 zlower = -1 * fmax(zih, zfh);
2910 zupper = -1 * fmin(zih, zfh);
2912 for (
int j = 0;
j < nMapComponents;
j++)
2914 hSeparatedMapComponent[
i][
j] =
new TH3F(Form(
"hIntDistortion%s_%s", sepAxis[
j].
Data(), side[
i].
Data()),
2915 Form(
"Integrated %s Deflection drifting from (phi,r,z) to z=endcap;phi;r;z (%s side)", sepAxis[
j].
Data(), side[
i].
Data()),
2916 nph, pih, pfh, nrh, rih, rfh, nzh, zlower, zupper);
2918 hValidToStepComponent[
i] =
new TH3F(Form(
"hValidToStep_%s", side[
i].
Data()),
2919 Form(
"Step before out of bounds drifting from (phi,r,z) to z=endcap;phi;r;z (%s side)", side[
i].
Data()),
2920 nph, pih, pfh, nrh, rih, rfh, nzh, zlower, zupper);
2921 hReachesReadout[
i] =
new TH3F(Form(
"hReachesReadout_%s", side[
i].
Data()),
2922 Form(
"Bool value for checking if successfuly drifting from (phi,r,z) to z=endcap and reaches the readout;phi;r;z (%s side)", side[
i].
Data()),
2923 nph, pih, pfh, nrh, rih, rfh, nzh, zlower, zupper);
2927 hRdeltaRComponent =
new TH2F(
"hRdeltaR",
"Bool value of input for a 2D R_versus_deltaR plot from Zstart to Zend;R;deltaR", nrh, rih, rfh, 10 * nrh, rih - 70, rfh + 170);
2934 TVector3
pos(((
int) (nrh / 2) + 0.5) * steplocal.Perp() + rih, 0,
zmin + ((int) (
nz / 2) + 0.5) *
step.Z());
2935 float posphi = ((int) (nph / 2) + 0.5) * steplocal.Phi() + pih;
2938 int xi[3] = {(int) floor((
pos.Perp() - rih) / steplocal.Perp()), (
int) floor((posphi - pih) / steplocal.Phi()), (
int) floor((
pos.Z() - zih) / steplocal.Z())};
2941 printf(
"rpz slice indices= (%d,%d,%d) (no twin)\n", xi[0], xi[1], xi[2]);
2943 int twinz = (-
pos.Z() - zih) / steplocal.Z();
2946 printf(
"rpz slice indices= (%d,%d,%d) twinz=%d\n", xi[0], xi[1], xi[2], twinz);
2949 const char axname[] =
"rpzrpz";
2950 int axn[] = {nrh, nph, nzh, nrh, nph, nzh};
2951 float axval[] = {(float)
pos.Perp(), (float)
pos.Phi(), (float)
pos.Z(), (float)
pos.Perp(), (float)
pos.Phi(), (float)
pos.Z()};
2952 float axbot[] = {rih, pih, zih, rih, pih, zih};
2953 float axtop[] = {rfh, pfh, zfh, rfh, pfh, zfh};
2954 TH2F *hIntDist[3][3];
2956 TH2F *hDiffDist[3][3];
2957 TH1F *hRDiffDist[2][3];
2958 for (
int i = 0;
i < 3;
i++)
2961 for (
int ax = 0;
ax < 3;
ax++)
2964 hDiffDist[
ax][
i] =
new TH2F(TString::Format(
"hDiffDist%c_%c%c", axname[
i], axname[
ax + 1], axname[
ax + 2]),
2965 TString::Format(
"%c component of diff. distortion in %c%c plane at %c=%2.3f;%c;%c",
2966 axname[i], axname[
ax + 1], axname[
ax + 2], axname[
ax], axval[ax], axname[ax + 1], axname[ax + 2]),
2967 axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2968 axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2969 hIntDist[
ax][
i] =
new TH2F(TString::Format(
"hIntDist%c_%c%c", axname[i], axname[ax + 1], axname[ax + 2]),
2970 TString::Format(
"%c component of int. distortion in %c%c plane at %c=%2.3f;%c;%c",
2971 axname[i], axname[ax + 1], axname[ax + 2], axname[ax], axval[ax], axname[ax + 1], axname[ax + 2]),
2972 axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
2973 axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
2975 hRDist[0][
i] =
new TH1F(TString::Format(
"hRDist%c", axname[
i]),
2976 TString::Format(
"%c component of int. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
2977 axname[i], axname[1], axval[1], axname[2], axval[2]),
2978 axn[0], axbot[0], axtop[0]);
2979 hRDist[1][
i] =
new TH1F(TString::Format(
"hRDist%cNeg", axname[i]),
2980 TString::Format(
"%c component of int. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
2981 axname[i], axname[1], axval[1], axname[2], axval[2]),
2982 axn[0], axbot[0], axtop[0]);
2983 hRDiffDist[0][
i] =
new TH1F(TString::Format(
"hRDiffDist%c", axname[i]),
2984 TString::Format(
"%c component of diff. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
2985 axname[i], axname[1], axval[1], axname[2], axval[2]),
2986 axn[0], axbot[0], axtop[0]);
2987 hRDiffDist[1][
i] =
new TH1F(TString::Format(
"hRDiffDist%cNeg", axname[i]),
2988 TString::Format(
"%c component of diff. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
2989 axname[i], axname[1], axval[1], axname[2], axval[2]),
2990 axn[0], axbot[0], axtop[0]);
2993 TVector3 inpart, outpart;
2994 TVector3 diffdistort, distort;
2999 float partR, partP, partZ;
3001 float distortR, distortP, distortZ;
3002 float distortX, distortY;
3003 float diffdistR, diffdistP, diffdistZ;
3004 TTree *dTree =
new TTree(
"dTree",
"Distortion per step z");
3005 dTree->Branch(
"r", &partR);
3006 dTree->Branch(
"p", &partP);
3007 dTree->Branch(
"z", &partZ);
3008 dTree->Branch(
"ir", &ir);
3009 dTree->Branch(
"ip", &ip);
3010 dTree->Branch(
"iz", &iz);
3011 dTree->Branch(
"dr", &distortR);
3012 dTree->Branch(
"drp", &distortP);
3013 dTree->Branch(
"dz", &distortZ);
3015 printf(
"generating separated distortion map with (%dx%dx%d) grid \n", nrh, nph, nzh);
3016 unsigned long long totalelements = nrh;
3017 totalelements *= nph;
3018 totalelements *= nzh;
3024 unsigned long long percent = totalelements / 100;
3026 printf(
"total elements = %llu\n", totalelements);
3037 inpart.SetXYZ(1, 0, 0);
3038 for (ir = 0; ir < nrh; ir++)
3040 partR = (ir + 0.5) * deltar + rih;
3043 inpart.SetPerp(partR + deltar);
3045 else if (ir == nrh - 1)
3047 inpart.SetPerp(partR - deltar);
3051 inpart.SetPerp(partR);
3053 for (ip = 0; ip < nph; ip++)
3055 partP = (ip + 0.5) * deltap + pih;
3056 inpart.SetPhi(partP);
3058 for (iz = 0; iz < nzh; iz++)
3060 partZ = (iz) *deltaz + zih;
3063 inpart.SetZ(partZ + deltaz);
3065 else if (iz == nzh - 1)
3067 inpart.SetZ(partZ - deltaz);
3073 partZ += 0.5 * deltaz;
3074 for (
int localside = 0; localside < nSides; localside++)
3079 distort =
GetTotalDistortion(z_readout, inpart, nSteps,
true, &validToStep, &successCheck);
3086 inpart.SetZ(-1 * inpart.Z());
3091 diffdistort.RotateZ(-inpart.Phi());
3092 diffdistP = diffdistort.Y();
3093 diffdistR = diffdistort.X();
3094 diffdistZ = diffdistort.Z();
3096 distortX = distort.X();
3097 distortY = distort.Y();
3098 distort.RotateZ(-inpart.Phi());
3099 distortP = distort.Y();
3100 distortR = distort.X();
3101 distortZ = distort.Z();
3103 float distComp[nMapComponents];
3104 distComp[0] = distortX;
3105 distComp[1] = distortY;
3106 distComp[2] = distortZ;
3107 distComp[3] = distortR;
3108 distComp[4] = distortP / partR;
3109 distComp[5] = distortP;
3111 for (
int c = 0;
c < nMapComponents;
c++)
3113 hSeparatedMapComponent[localside][
c]->Fill(partP, partR, partZ, distComp[
c]);
3116 hValidToStepComponent[localside]->Fill(partP, partR, partZ, validToStep);
3117 hReachesReadout[localside]->Fill(partP, partR, partZ, successCheck);
3123 hIntDistortionR->Fill(partP, partR, partZ, distortR);
3124 hIntDistortionP->Fill(partP, partR, partZ, distortP);
3125 hIntDistortionZ->Fill(partP, partR, partZ, distortZ);
3129 hIntDistortionX->Fill(partP, partR, partZ, distortX);
3130 hIntDistortionY->Fill(partP, partR, partZ, distortY);
3134 if (ir == xi[0] && localside == 0)
3137 hIntDist[0][0]->Fill(partP, partZ, distortR);
3138 hIntDist[0][1]->Fill(partP, partZ, distortP);
3139 hIntDist[0][2]->Fill(partP, partZ, distortZ);
3140 hDiffDist[0][0]->Fill(partP, partZ, diffdistR);
3141 hDiffDist[0][1]->Fill(partP, partZ, diffdistP);
3142 hDiffDist[0][2]->Fill(partP, partZ, diffdistZ);
3144 if (ip == xi[1] && localside == 0)
3147 hIntDist[1][0]->Fill(partZ, partR, distortR);
3148 hIntDist[1][1]->Fill(partZ, partR, distortP);
3149 hIntDist[1][2]->Fill(partZ, partR, distortZ);
3150 hDiffDist[1][0]->Fill(partZ, partR, diffdistR);
3151 hDiffDist[1][1]->Fill(partZ, partR, diffdistP);
3152 hDiffDist[1][2]->Fill(partZ, partR, diffdistZ);
3154 if (iz == xi[2] && localside == 0)
3156 hRDist[0][0]->Fill(partR, distortR);
3157 hRDist[0][1]->Fill(partR, distortP);
3158 hRDist[0][2]->Fill(partR, distortZ);
3159 hRDiffDist[0][0]->Fill(partR, diffdistR);
3160 hRDiffDist[0][1]->Fill(partR, diffdistP);
3161 hRDiffDist[0][2]->Fill(partR, diffdistZ);
3163 if (
hasTwin && iz == twinz && localside == 1)
3165 hRDist[1][0]->Fill(partR, distortR);
3166 hRDist[1][1]->Fill(partR, distortP);
3167 hRDist[1][2]->Fill(partR, distortZ);
3168 hRDiffDist[1][0]->Fill(partR, diffdistR);
3169 hRDiffDist[1][1]->Fill(partR, diffdistP);
3170 hRDiffDist[1][2]->Fill(partR, diffdistZ);
3173 if (iz == xi[2] && localside == 0)
3177 hIntDist[2][0]->Fill(partR, partP, distortR);
3178 hIntDist[2][1]->Fill(partR, partP, distortP);
3179 hIntDist[2][2]->Fill(partR, partP, distortZ);
3180 hDiffDist[2][0]->Fill(partR, partP, diffdistR);
3181 hDiffDist[2][1]->Fill(partR, partP, diffdistP);
3182 hDiffDist[2][2]->Fill(partR, partP, diffdistZ);
3185 if (!(el % waypoint))
3187 printf(
"generating distortions %d%%: ", (
int) (el / percent));
3188 printf(
"distortion at (ir=%d,ip=%d,iz=%d) is (%E,%E,%E)\n",
3189 ir, ip, iz, distortR, distortP, distortZ);
3196 printf(
"Completed distortion generation. Saving outputs...\n");
3198 TCanvas *canvas =
new TCanvas(
"cdistort",
"distortion integrals", 1200, 800);
3200 printf(
"was able to make a tcanvas\n");
3202 TPad *
c =
new TPad(
"cplots",
"distortion integral plots", 0, 0.2, 1, 1);
3204 TPad *textpad =
new TPad(
"ctext",
"distortion integral plots", 0, 0.0, 1, 0.2);
3205 printf(
"was able to make some tpads\n");
3208 gStyle->SetOptStat();
3209 printf(
"was able to interact with gStyle\n");
3211 for (
int i = 0;
i < 3;
i++)
3214 for (
int ax = 0;
ax < 3;
ax++)
3216 printf(
"looping over components i=%d ax=%d\n",
i,
ax);
3219 c->cd(
i * 4 +
ax + 1);
3220 gPad->SetRightMargin(0.15);
3221 hIntDist[
ax][
i]->SetStats(
false);
3222 hIntDist[
ax][
i]->Draw(
"colz");
3225 printf(
"drawing R profile %d\n",
i);
3228 hRDist[0][
i]->SetStats(
false);
3229 hRDist[0][
i]->SetFillColor(kRed);
3230 hRDist[0][
i]->Draw(
"hist");
3233 printf(
"drawing R profile twin %d\n",
i);
3234 hRDist[1][
i]->SetStats(
false);
3235 hRDist[1][
i]->SetLineColor(kBlue);
3236 hRDist[1][
i]->Draw(
"hist,same");
3239 printf(
"switching to textpad\n");
3243 float texshift = 0.12;
3244 TLatex *
tex =
new TLatex(0.0, texpos,
"Fill Me In");
3245 printf(
"built TLatex\n");
3247 tex->SetTextSize(texshift * 0.8);
3253 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3265 printf(
"cd'ing to canvas:\n");
3274 printf(
"was able to complete drawing on both pads\n");
3276 canvas->SaveAs(summaryFilename.Data());
3284 for (
int i = 0;
i < 3;
i++)
3287 for (
int ax = 0;
ax < 3;
ax++)
3290 c->cd(
i * 4 +
ax + 1);
3291 gPad->SetRightMargin(0.15);
3292 hDiffDist[
ax][
i]->SetStats(
false);
3293 hDiffDist[
ax][
i]->Draw(
"colz");
3296 hRDiffDist[0][
i]->SetStats(
false);
3297 hRDiffDist[0][
i]->SetFillColor(kRed);
3298 hRDiffDist[0][
i]->Draw(
"hist");
3301 hRDiffDist[1][
i]->SetStats(
false);
3302 hRDiffDist[1][
i]->SetLineColor(kBlue);
3303 hRDiffDist[1][
i]->Draw(
"hist same");
3309 tex->SetTextSize(texshift * 0.8);
3315 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3321 tex->DrawLatex(0.05, texpos,
"Differential Plots");
3334 canvas->SaveAs(diffSummaryFilename.Data());
3336 printf(
"saving map histograms to:%s.\n", distortionFilename.Data());
3339 for (
int i = 0;
i < nSides;
i++)
3342 for (
int j = 0;
j < nMapComponents;
j++)
3344 hSeparatedMapComponent[
i][
j]->GetSumw2()->Set(0);
3345 hSeparatedMapComponent[
i][
j]->Write();
3347 hValidToStepComponent[
i]->GetSumw2()->Set(0);
3348 hValidToStepComponent[
i]->Write();
3350 hReachesReadout[
i]->GetSumw2()->Set(0);
3351 hReachesReadout[
i]->Write();
3359 hDistortionR->GetSumw2()->Set(0);
3360 hDistortionP->GetSumw2()->Set(0);
3361 hDistortionZ->GetSumw2()->Set(0);
3362 hIntDistortionR->GetSumw2()->Set(0);
3363 hIntDistortionP->GetSumw2()->Set(0);
3364 hIntDistortionZ->GetSumw2()->Set(0);
3367 hIntDistortionX->GetSumw2()->Set(0);
3368 hIntDistortionY->GetSumw2()->Set(0);
3370 hDistortionR->Write();
3371 hDistortionP->Write();
3372 hDistortionZ->Write();
3373 hIntDistortionR->Write();
3374 hIntDistortionP->Write();
3375 hIntDistortionZ->Write();
3376 if (
false && andCartesian)
3378 hIntDistortionX->Write();
3379 hIntDistortionY->Write();
3381 printf(
"finished writing histograms\n");
3387 printf(
"wrote separated map and summary to %s.\n", filebase);
3388 printf(
"map:%s.\n", distortionFilename.Data());
3389 printf(
"summary:%s.\n", summaryFilename.Data());
3399 printf(
"WARNING: You called the version of the distortion generator that generates a unified map. Are you sure you meant to do this?\n");
3402 bool makeUnifiedMap =
true;
3404 TVector3 steplocal(
step.Perp() / r_subsamples, 0,
step.Z() / z_subsamples);
3405 steplocal.SetPhi(
step.Phi() / p_subsamples);
3406 float deltar = steplocal.Perp();
3407 float deltap = steplocal.Phi();
3408 float deltaz = steplocal.Z();
3409 TVector3 stepzvec(0, 0, deltaz);
3423 int nph =
nphi * p_subsamples + 2;
3424 int nrh =
nr * r_subsamples + 2;
3425 int nzh =
nz * z_subsamples + 2;
3428 if (
hasTwin && makeUnifiedMap)
3430 lowerEdge.SetZ(-1 * upperEdge.Z());
3431 nzh +=
nz * z_subsamples;
3434 float rih = lowerEdge.Perp() - 0.5 *
step.Perp() - steplocal.Perp();
3435 float rfh = upperEdge.Perp() + 0.5 *
step.Perp() + steplocal.Perp();
3436 float pih =
FilterPhiPos(lowerEdge.Phi()) - 0.5 *
step.Phi() - steplocal.Phi();
3437 float pfh =
FilterPhiPos(upperEdge.Phi()) + 0.5 *
step.Phi() + steplocal.Phi();
3438 float zih = lowerEdge.Z() - 0.5 *
step.Z() - steplocal.Z();
3439 float zfh = upperEdge.Z() + 0.5 *
step.Z() + steplocal.Z();
3440 float z_readout = upperEdge.Z() + 0.5 *
step.Z();
3442 printf(
"generating distortion map...\n");
3443 printf(
"file=%s\n", filebase);
3447 TString distortionFilename;
3448 distortionFilename.Form(
"%s.distortion_map.hist.root", filebase);
3449 TString summaryFilename;
3450 summaryFilename.Form(
"%s.distortion_summary.pdf", filebase);
3451 TString diffSummaryFilename;
3452 diffSummaryFilename.Form(
"%s.differential_summary.pdf", filebase);
3454 TFile *outf = TFile::Open(distortionFilename.Data(),
"RECREATE");
3459 TH3F *hDistortionR =
new TH3F(
"hDistortionR",
"Per-z-bin Distortion in the R direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3460 TH3F *hDistortionP =
new TH3F(
"hDistortionP",
"Per-z-bin Distortion in the RPhi direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3461 TH3F *hDistortionZ =
new TH3F(
"hDistortionZ",
"Per-z-bin Distortion in the Z direction as a function of (phi,r,z) (centered in r,phi, z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3462 TH3F *hIntDistortionR =
new TH3F(
"hIntDistortionR",
"Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3463 TH3F *hIntDistortionP =
new TH3F(
"hIntDistortionP",
"Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3464 TH3F *hIntDistortionZ =
new TH3F(
"hIntDistortionZ",
"Integrated R Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3466 TH3F *hIntDistortionX =
new TH3F(
"hIntDistortionX",
"Integrated X Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3467 TH3F *hIntDistortionY =
new TH3F(
"hIntDistortionY",
"Integrated Y Distortion from (phi,r,z) to z=0 (centered in r,phi, and z);phi;r;z", nph, pih, pfh, nrh, rih, rfh, nzh, zih, zfh);
3482 TVector3
pos(((
int) (nrh / 2) + 0.5) * steplocal.Perp() + rih, 0,
zmin + ((int) (
nz / 2) + 0.5) *
step.Z());
3483 float posphi = ((int) (nph / 2) + 0.5) * steplocal.Phi() + pih;
3486 int xi[3] = {(int) floor((
pos.Perp() - rih) / steplocal.Perp()), (
int) floor((posphi - pih) / steplocal.Phi()), (
int) floor((
pos.Z() - zih) / steplocal.Z())};
3489 printf(
"rpz slice indices= (%d,%d,%d) (no twin)\n", xi[0], xi[1], xi[2]);
3495 printf(
"rpz slice indices= (%d,%d,%d) twinz=%d\n", xi[0], xi[1], xi[2], twinz);
3498 const char axname[] =
"rpzrpz";
3499 int axn[] = {nrh, nph, nzh, nrh, nph, nzh};
3500 float axval[] = {(float)
pos.Perp(), (float)
pos.Phi(), (float)
pos.Z(), (float)
pos.Perp(), (float)
pos.Phi(), (float)
pos.Z()};
3501 float axbot[] = {rih, pih, zih, rih, pih, zih};
3502 float axtop[] = {rfh, pfh, zfh, rfh, pfh, zfh};
3503 TH2F *hIntDist[3][3];
3505 TH2F *hDiffDist[3][3];
3506 TH1F *hRDiffDist[2][3];
3507 for (
int i = 0;
i < 3;
i++)
3510 for (
int ax = 0;
ax < 3;
ax++)
3513 hDiffDist[
ax][
i] =
new TH2F(TString::Format(
"hDiffDist%c_%c%c", axname[
i], axname[
ax + 1], axname[
ax + 2]),
3514 TString::Format(
"%c component of diff. distortion in %c%c plane at %c=%2.3f;%c;%c",
3515 axname[i], axname[
ax + 1], axname[
ax + 2], axname[
ax], axval[ax], axname[ax + 1], axname[ax + 2]),
3516 axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
3517 axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
3518 hIntDist[
ax][
i] =
new TH2F(TString::Format(
"hIntDist%c_%c%c", axname[i], axname[ax + 1], axname[ax + 2]),
3519 TString::Format(
"%c component of int. distortion in %c%c plane at %c=%2.3f;%c;%c",
3520 axname[i], axname[ax + 1], axname[ax + 2], axname[ax], axval[ax], axname[ax + 1], axname[ax + 2]),
3521 axn[ax + 1], axbot[ax + 1], axtop[ax + 1],
3522 axn[ax + 2], axbot[ax + 2], axtop[ax + 2]);
3524 hRDist[0][
i] =
new TH1F(TString::Format(
"hRDist%c", axname[
i]),
3525 TString::Format(
"%c component of int. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
3526 axname[i], axname[1], axval[1], axname[2], axval[2]),
3527 axn[0], axbot[0], axtop[0]);
3528 hRDist[1][
i] =
new TH1F(TString::Format(
"hRDist%cNeg", axname[i]),
3529 TString::Format(
"%c component of int. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
3530 axname[i], axname[1], axval[1], axname[2], axval[2]),
3531 axn[0], axbot[0], axtop[0]);
3532 hRDiffDist[0][
i] =
new TH1F(TString::Format(
"hRDist%c", axname[i]),
3533 TString::Format(
"%c component of diff. distortion vs r with %c=%2.3f and %c=%2.3f;r(cm);#delta (cm)",
3534 axname[i], axname[1], axval[1], axname[2], axval[2]),
3535 axn[0], axbot[0], axtop[0]);
3536 hRDiffDist[1][
i] =
new TH1F(TString::Format(
"hRDist%cNeg", axname[i]),
3537 TString::Format(
"%c component of diff. distortion vs r with %c=%2.3f and %c=-%2.3f;r(cm);#delta (cm)",
3538 axname[i], axname[1], axval[1], axname[2], axval[2]),
3539 axn[0], axbot[0], axtop[0]);
3542 TVector3 inpart, outpart;
3547 float partR, partP, partZ;
3549 float distortR, distortP, distortZ;
3550 float distortX, distortY;
3551 float diffdistR, diffdistP, diffdistZ;
3552 TTree *dTree =
new TTree(
"dTree",
"Distortion per step z");
3553 dTree->Branch(
"r", &partR);
3554 dTree->Branch(
"p", &partP);
3555 dTree->Branch(
"z", &partZ);
3556 dTree->Branch(
"ir", &ir);
3557 dTree->Branch(
"ip", &ip);
3558 dTree->Branch(
"iz", &iz);
3559 dTree->Branch(
"dr", &distortR);
3560 dTree->Branch(
"drp", &distortP);
3561 dTree->Branch(
"dz", &distortZ);
3563 printf(
"generating distortion map with (%dx%dx%d) grid \n", nrh, nph, nzh);
3564 unsigned long long totalelements = nrh;
3565 totalelements *= nph;
3566 totalelements *= nzh;
3567 unsigned long long percent = totalelements / 100 *
debug_npercent;
3568 printf(
"total elements = %llu\n", totalelements);
3579 inpart.SetXYZ(1, 0, 0);
3580 for (ir = 0; ir < nrh; ir++)
3582 partR = (ir + 0.5) * deltar + rih;
3585 inpart.SetPerp(partR + deltar);
3587 else if (ir == nrh - 1)
3589 inpart.SetPerp(partR - deltar);
3593 inpart.SetPerp(partR);
3595 for (ip = 0; ip < nph; ip++)
3597 partP = (ip + 0.5) * deltap + pih;
3598 inpart.SetPhi(partP);
3600 for (iz = 0; iz < nzh; iz++)
3602 partZ = (iz) *deltaz + zih;
3605 inpart.SetZ(partZ + deltaz);
3607 else if (iz == nzh - 1)
3609 inpart.SetZ(partZ - deltaz);
3615 partZ += 0.5 * deltaz;
3621 if (
hasTwin && inpart.Z() < 0)
3623 distort =
twin->
GetTotalDistortion(inpart.Z(), inpart + stepzvec, nSteps,
true, &validToStep, &successCheck);
3627 distort =
GetTotalDistortion(inpart.Z() + deltaz, inpart, nSteps,
true, &validToStep, &successCheck);
3629 distort.RotateZ(-inpart.Phi());
3630 diffdistP = distort.Y();
3631 diffdistR = distort.X();
3632 diffdistZ = distort.Z();
3633 hDistortionR->Fill(partP, partR, partZ, diffdistR);
3634 hDistortionP->Fill(partP, partR, partZ, diffdistP);
3635 hDistortionZ->Fill(partP, partR, partZ, diffdistZ);
3639 if (
hasTwin && makeUnifiedMap && inpart.Z() < 0)
3641 distort =
twin->
GetTotalDistortion(-z_readout, inpart + stepzvec, nSteps,
true, &validToStep, &successCheck);
3645 distort =
GetTotalDistortion(z_readout, inpart, nSteps,
true, &validToStep, &successCheck);
3647 distortX = distort.X();
3648 distortY = distort.Y();
3649 distort.RotateZ(-inpart.Phi());
3650 distortP = distort.Y();
3651 distortR = distort.X();
3652 distortZ = distort.Z();
3658 hIntDistortionR->Fill(partP, partR, partZ, distortR);
3659 hIntDistortionP->Fill(partP, partR, partZ, distortP);
3660 hIntDistortionZ->Fill(partP, partR, partZ, distortZ);
3664 hIntDistortionX->Fill(partP, partR, partZ, distortX);
3665 hIntDistortionY->Fill(partP, partR, partZ, distortY);
3672 hIntDist[0][0]->Fill(partP, partZ, distortR);
3673 hIntDist[0][1]->Fill(partP, partZ, distortP);
3674 hIntDist[0][2]->Fill(partP, partZ, distortZ);
3675 hDiffDist[0][0]->Fill(partP, partZ, diffdistR);
3676 hDiffDist[0][1]->Fill(partP, partZ, diffdistP);
3677 hDiffDist[0][2]->Fill(partP, partZ, diffdistZ);
3682 hIntDist[1][0]->Fill(partZ, partR, distortR);
3683 hIntDist[1][1]->Fill(partZ, partR, distortP);
3684 hIntDist[1][2]->Fill(partZ, partR, distortZ);
3685 hDiffDist[1][0]->Fill(partZ, partR, diffdistR);
3686 hDiffDist[1][1]->Fill(partZ, partR, diffdistP);
3687 hDiffDist[1][2]->Fill(partZ, partR, diffdistZ);
3691 hRDist[0][0]->Fill(partR, distortR);
3692 hRDist[0][1]->Fill(partR, distortP);
3693 hRDist[0][2]->Fill(partR, distortZ);
3694 hRDiffDist[0][0]->Fill(partR, diffdistR);
3695 hRDiffDist[0][1]->Fill(partR, diffdistP);
3696 hRDiffDist[0][2]->Fill(partR, diffdistZ);
3700 hRDist[1][0]->Fill(partR, distortR);
3701 hRDist[1][1]->Fill(partR, distortP);
3702 hRDist[1][2]->Fill(partR, distortZ);
3703 hRDiffDist[1][0]->Fill(partR, diffdistR);
3704 hRDiffDist[1][1]->Fill(partR, diffdistP);
3705 hRDiffDist[1][2]->Fill(partR, diffdistZ);
3712 hIntDist[2][0]->Fill(partR, partP, distortR);
3713 hIntDist[2][1]->Fill(partR, partP, distortP);
3714 hIntDist[2][2]->Fill(partR, partP, distortZ);
3715 hDiffDist[2][0]->Fill(partR, partP, diffdistR);
3716 hDiffDist[2][1]->Fill(partR, partP, diffdistP);
3717 hDiffDist[2][2]->Fill(partR, partP, diffdistZ);
3720 if (!(el % percent))
3723 printf(
"distortion at (ir=%d,ip=%d,iz=%d) is (%E,%E,%E)\n",
3724 ir, ip, iz, distortR, distortP, distortZ);
3731 TCanvas *canvas =
new TCanvas(
"cdistort",
"distortion integrals", 1200, 800);
3734 TPad *
c =
new TPad(
"cplots",
"distortion integral plots", 0, 0.2, 1, 1);
3736 TPad *textpad =
new TPad(
"ctext",
"distortion integral plots", 0, 0.0, 1, 0.2);
3738 gStyle->SetOptStat();
3739 for (
int i = 0;
i < 3;
i++)
3742 for (
int ax = 0;
ax < 3;
ax++)
3745 c->cd(
i * 4 +
ax + 1);
3746 gPad->SetRightMargin(0.15);
3747 hIntDist[
ax][
i]->SetStats(
false);
3748 hIntDist[
ax][
i]->Draw(
"colz");
3751 hRDist[0][
i]->SetStats(
false);
3752 hRDist[0][
i]->SetFillColor(kRed);
3753 hRDist[0][
i]->Draw(
"hist");
3756 hRDist[1][
i]->SetStats(
false);
3757 hRDist[1][
i]->SetLineColor(kBlue);
3758 hRDist[1][
i]->Draw(
"hist,same");
3777 float texshift = 0.12;
3778 TLatex *
tex =
new TLatex(0.0, texpos,
"Fill Me In");
3779 tex->SetTextSize(texshift * 0.8);
3785 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3802 canvas->SaveAs(summaryFilename.Data());
3810 for (
int i = 0;
i < 3;
i++)
3813 for (
int ax = 0;
ax < 3;
ax++)
3816 c->cd(
i * 4 +
ax + 1);
3817 gPad->SetRightMargin(0.15);
3818 hDiffDist[
ax][
i]->SetStats(
false);
3819 hDiffDist[
ax][
i]->Draw(
"colz");
3822 hRDiffDist[0][
i]->SetStats(
false);
3823 hRDiffDist[0][
i]->SetFillColor(kRed);
3824 hRDiffDist[0][
i]->Draw(
"hist");
3827 hRDiffDist[1][
i]->SetStats(
false);
3828 hRDiffDist[1][
i]->SetLineColor(kBlue);
3829 hRDiffDist[1][
i]->Draw(
"hist same");
3835 tex->SetTextSize(texshift * 0.8);
3841 tex->DrawLatex(0.05, texpos, Form(
"Drifting grid of (rp)=(%d x %d) electrons with %d steps", nrh, nph, nSteps));
3847 tex->DrawLatex(0.05, texpos,
"Differential Plots");
3860 canvas->SaveAs(diffSummaryFilename.Data());
3866 hDistortionR->Write();
3867 hDistortionP->Write();
3868 hDistortionZ->Write();
3869 hIntDistortionR->Write();
3870 hIntDistortionP->Write();
3871 hIntDistortionZ->Write();
3874 hIntDistortionX->Write();
3875 hIntDistortionY->Write();
3902 printf(
"wrote map and summary to %s.\n", filebase);
3903 printf(
"map:%s.\n", distortionFilename.Data());
3904 printf(
"summary:%s.\n", summaryFilename.Data());
3911 int defaultsteps = 100;
3917 return swimToInSteps(zdest, start, defaultsteps, interpolate, &goodtostep);
3938 double zdist = zdest - start.Z();
3957 else if (interpolate)
3970 printf(
"GetStepDistortion is attempting to swim with no drift field:\n");
3971 printf(
"GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
3972 printf(
"GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
3976 double EfieldZ = fieldInt.Z() / zdist;
3977 double BfieldZ = fieldIntB.Z() / zdist;
3989 double c0 = 1 / (1 + T2om2);
3990 double c1 = T1om / (1 + T1om * T1om);
3991 double c2 = T2om2 / (1 + T2om2);
3993 TVector3 EintOverEz = 1 / EfieldZ * fieldInt;
3994 TVector3 BintOverBz = 1 / BfieldZ * fieldIntB;
3999 double deltaX = c0 * EintOverEz.X() + c1 * EintOverEz.Y() - c1 * BintOverBz.Y() + c2 * BintOverBz.X();
4000 double deltaY = c0 * EintOverEz.Y() - c1 * EintOverEz.X() + c1 * BintOverBz.X() + c2 * BintOverBz.Y();
4003 double vprime = (5000 *
cm /
s) / (100 *
V /
cm);
4004 double vdoubleprime = 0;
4009 double deltaZ = vprime /
vdrift * (fieldInt.Z() - zdist *
Enominal) + vdoubleprime /
vdrift * (fieldInt.Z() -
Enominal * zdist) * (fieldInt.Z() -
Enominal * zdist) / (2 * zdist) - 0.5 / zdist * (EintOverEz.X() * EintOverEz.X() + EintOverEz.Y() * EintOverEz.Y()) + c1 / zdist * (EintOverEz.X() * BintOverBz.Y() - EintOverEz.Y() * BintOverBz.X()) + c2 / zdist * (EintOverEz.X() * BintOverBz.X() + EintOverEz.Y() * BintOverBz.Y()) + c2 / zdist * (BintOverBz.X() * BintOverBz.X() + BintOverBz.Y() * BintOverBz.Y());
4013 printf(
"GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4014 printf(
"GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4015 printf(
"GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4016 printf(
"GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
4017 printf(
"GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4018 printf(
"GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4024 printf(
"GetStepDistortion produced a very large zdistortion!\n");
4025 printf(
"GetStepDistortion: zdist=%2.4f, deltaZ=%2.4f, Delta/z=%2.4f\n", zdist, deltaZ, deltaZ / zdist);
4026 printf(
"GetStepDistortion: average field Z: Ez=%2.4fV/cm, Bz=%2.4fT\n", EfieldZ / (
V /
cm), BfieldZ /
Tesla);
4027 printf(
"GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4028 printf(
"GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4029 printf(
"GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4030 printf(
"GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
4031 printf(
"GetStepDistortion: EfieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4032 printf(
"GetStepDistortion: BfieldInt=(%E,%E,%E)\n", fieldIntB.X(), fieldIntB.Y(), fieldIntB.Z());
4033 printf(
"GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4039 printf(
"GetStepDistortion produced a very small deltaX: %E\n", deltaX);
4040 printf(
"GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4041 printf(
"GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4042 printf(
"GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4043 printf(
"GetStepDistortion: (%2.4f,%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), zdest);
4044 printf(
"GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4045 printf(
"GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4052 printf(
"GetStepDistortion produced a very large deltaX: %E\n", deltaX);
4053 printf(
"GetStepDistortion: (c0,c1,c2)=(%E,%E,%E)\n", c0, c1, c2);
4054 printf(
"GetStepDistortion: EintOverEz==(%E,%E,%E)\n", EintOverEz.X(), EintOverEz.Y(), EintOverEz.Z());
4055 printf(
"GetStepDistortion: BintOverBz==(%E,%E,%E)\n", BintOverBz.X(), BintOverBz.Y(), BintOverBz.Z());
4056 printf(
"GetStepDistortion: (%2.4f,%2.4f,%2.4f) (rp)=(%2.4f,%2.4f) to z=%2.4f\n", start.X(), start.Y(), start.Z(), start.Perp(), start.Phi(), zdest);
4057 printf(
"GetStepDistortion: fieldInt=(%E,%E,%E)\n", fieldInt.X(), fieldInt.Y(), fieldInt.Z());
4058 printf(
"GetStepDistortion: delta=(%E,%E,%E)\n", deltaX, deltaY, deltaZ);
4064 TVector3 shift(deltaX, deltaY, deltaZ);
4067 shift.RotateZ(-start.Phi());
4071 shift.RotateZ(start.Phi());
4162 printf(
"Caution: tried to read charge at zbin=%d! No twin available to handle this\n", z);