LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
SpaceChargeStandard.cxx
Go to the documentation of this file.
1 // \file SpaceChargeStandard.cxx
3 //
4 // \brief implementation of class for storing/accessing space charge distortions (exists as an example, by default this does nothing, instead requiring experiment-specific implementation)
5 //
6 // \author mrmooney@bnl.gov
7 //
9 
10 // C++ language includes
11 #include <iostream>
12 #include <fstream>
13 #include <string>
14 #include <vector>
15 #include "math.h"
16 #include "stdio.h"
17 
18 // LArSoft includes
20 
21 // Framework includes
23 
24 //-----------------------------------------------
26  fhicl::ParameterSet const& pset
27 )
28 {
29  Configure(pset);
30 }
31 
32 //------------------------------------------------
34 {
35  fEnableSimSpatialSCE = pset.get<bool>("EnableSimSpatialSCE");
36  fEnableSimEfieldSCE = pset.get<bool>("EnableSimEfieldSCE");
37  fEnableCorrSCE = pset.get<bool>("EnableCorrSCE");
38 
39  // check that the old obsoleted parameter is not in configuration:
40  if (pset.has_key("EnableSimulationSCE")) {
42  << "Configuration parameter 'EnableSimulationSCE' has been replaced by 'EnableSimSpatialSCE'.\n";
43  }
44 
45  if((fEnableSimSpatialSCE == true) | (fEnableSimEfieldSCE == true))
46  {
47  fRepresentationType = pset.get<std::string>("RepresentationType");
48  fInputFilename = pset.get<std::string>("InputFilename");
49 
50  std::string fname;
51  cet::search_path sp("FW_SEARCH_PATH");
52  sp.find_file(fInputFilename,fname);
53 
54  std::unique_ptr<TFile> infile(new TFile(fname.c_str(), "READ"));
55  if(!infile->IsOpen()) throw art::Exception(art::errors::Configuration) << "Could not find the space charge effect file '" << fname << "'!\n";
56 
57  if(fRepresentationType == "Parametric")
58  {
59  for(int i = 0; i < 5; i++)
60  {
61  g1_x[i] = (TGraph*)infile->Get(Form("deltaX/g1_%d",i));
62  g2_x[i] = (TGraph*)infile->Get(Form("deltaX/g2_%d",i));
63  g3_x[i] = (TGraph*)infile->Get(Form("deltaX/g3_%d",i));
64  g4_x[i] = (TGraph*)infile->Get(Form("deltaX/g4_%d",i));
65  g5_x[i] = (TGraph*)infile->Get(Form("deltaX/g5_%d",i));
66 
67  g1_y[i] = (TGraph*)infile->Get(Form("deltaY/g1_%d",i));
68  g2_y[i] = (TGraph*)infile->Get(Form("deltaY/g2_%d",i));
69  g3_y[i] = (TGraph*)infile->Get(Form("deltaY/g3_%d",i));
70  g4_y[i] = (TGraph*)infile->Get(Form("deltaY/g4_%d",i));
71  g5_y[i] = (TGraph*)infile->Get(Form("deltaY/g5_%d",i));
72  g6_y[i] = (TGraph*)infile->Get(Form("deltaY/g6_%d",i));
73 
74  g1_z[i] = (TGraph*)infile->Get(Form("deltaZ/g1_%d",i));
75  g2_z[i] = (TGraph*)infile->Get(Form("deltaZ/g2_%d",i));
76  g3_z[i] = (TGraph*)infile->Get(Form("deltaZ/g3_%d",i));
77  g4_z[i] = (TGraph*)infile->Get(Form("deltaZ/g4_%d",i));
78 
79  g1_Ex[i] = (TGraph*)infile->Get(Form("deltaExOverE/g1_%d",i));
80  g2_Ex[i] = (TGraph*)infile->Get(Form("deltaExOverE/g2_%d",i));
81  g3_Ex[i] = (TGraph*)infile->Get(Form("deltaExOverE/g3_%d",i));
82  g4_Ex[i] = (TGraph*)infile->Get(Form("deltaExOverE/g4_%d",i));
83  g5_Ex[i] = (TGraph*)infile->Get(Form("deltaExOverE/g5_%d",i));
84 
85  g1_Ey[i] = (TGraph*)infile->Get(Form("deltaEyOverE/g1_%d",i));
86  g2_Ey[i] = (TGraph*)infile->Get(Form("deltaEyOverE/g2_%d",i));
87  g3_Ey[i] = (TGraph*)infile->Get(Form("deltaEyOverE/g3_%d",i));
88  g4_Ey[i] = (TGraph*)infile->Get(Form("deltaEyOverE/g4_%d",i));
89  g5_Ey[i] = (TGraph*)infile->Get(Form("deltaEyOverE/g5_%d",i));
90  g6_Ey[i] = (TGraph*)infile->Get(Form("deltaEyOverE/g6_%d",i));
91 
92  g1_Ez[i] = (TGraph*)infile->Get(Form("deltaEzOverE/g1_%d",i));
93  g2_Ez[i] = (TGraph*)infile->Get(Form("deltaEzOverE/g2_%d",i));
94  g3_Ez[i] = (TGraph*)infile->Get(Form("deltaEzOverE/g3_%d",i));
95  g4_Ez[i] = (TGraph*)infile->Get(Form("deltaEzOverE/g4_%d",i));
96  }
97 
98  g1_x[5] = (TGraph*)infile->Get("deltaX/g1_5");
99  g2_x[5] = (TGraph*)infile->Get("deltaX/g2_5");
100  g3_x[5] = (TGraph*)infile->Get("deltaX/g3_5");
101  g4_x[5] = (TGraph*)infile->Get("deltaX/g4_5");
102  g5_x[5] = (TGraph*)infile->Get("deltaX/g5_5");
103 
104  g1_y[5] = (TGraph*)infile->Get("deltaY/g1_5");
105  g2_y[5] = (TGraph*)infile->Get("deltaY/g2_5");
106  g3_y[5] = (TGraph*)infile->Get("deltaY/g3_5");
107  g4_y[5] = (TGraph*)infile->Get("deltaY/g4_5");
108  g5_y[5] = (TGraph*)infile->Get("deltaY/g5_5");
109  g6_y[5] = (TGraph*)infile->Get("deltaY/g6_5");
110 
111  g1_x[6] = (TGraph*)infile->Get("deltaX/g1_6");
112  g2_x[6] = (TGraph*)infile->Get("deltaX/g2_6");
113  g3_x[6] = (TGraph*)infile->Get("deltaX/g3_6");
114  g4_x[6] = (TGraph*)infile->Get("deltaX/g4_6");
115  g5_x[6] = (TGraph*)infile->Get("deltaX/g5_6");
116 
117  g1_Ex[5] = (TGraph*)infile->Get("deltaExOverE/g1_5");
118  g2_Ex[5] = (TGraph*)infile->Get("deltaExOverE/g2_5");
119  g3_Ex[5] = (TGraph*)infile->Get("deltaExOverE/g3_5");
120  g4_Ex[5] = (TGraph*)infile->Get("deltaExOverE/g4_5");
121  g5_Ex[5] = (TGraph*)infile->Get("deltaExOverE/g5_5");
122 
123  g1_Ey[5] = (TGraph*)infile->Get("deltaEyOverE/g1_5");
124  g2_Ey[5] = (TGraph*)infile->Get("deltaEyOverE/g2_5");
125  g3_Ey[5] = (TGraph*)infile->Get("deltaEyOverE/g3_5");
126  g4_Ey[5] = (TGraph*)infile->Get("deltaEyOverE/g4_5");
127  g5_Ey[5] = (TGraph*)infile->Get("deltaEyOverE/g5_5");
128  g6_Ey[5] = (TGraph*)infile->Get("deltaEyOverE/g6_5");
129 
130  g1_Ex[6] = (TGraph*)infile->Get("deltaExOverE/g1_6");
131  g2_Ex[6] = (TGraph*)infile->Get("deltaExOverE/g2_6");
132  g3_Ex[6] = (TGraph*)infile->Get("deltaExOverE/g3_6");
133  g4_Ex[6] = (TGraph*)infile->Get("deltaExOverE/g4_6");
134  g5_Ex[6] = (TGraph*)infile->Get("deltaExOverE/g5_6");
135  }
136 
137  infile->Close();
138  }
139 
140  if(fEnableCorrSCE == true)
141  {
142  // Grab other parameters from pset
143  }
144 
145  return true;
146 }
147 
148 //------------------------------------------------
150 {
151  if (ts == 0) return false;
152 
153  return true;
154 }
155 
156 //----------------------------------------------------------------------------
160 {
161  return fEnableSimSpatialSCE;
162 }
163 
164 //----------------------------------------------------------------------------
168 {
169  return fEnableSimEfieldSCE;
170 }
171 
172 //----------------------------------------------------------------------------
175 {
176  return fEnableCorrSCE;
177 }
178 
179 //----------------------------------------------------------------------------
183 {
184  std::vector<double> thePosOffsets;
185 
186  if(IsInsideBoundaries(point.X(), point.Y(), point.Z()) == false)
187  {
188  thePosOffsets.resize(3,0.0);
189  }
190  else
191  {
192  if(fRepresentationType == "Parametric")
193  thePosOffsets = GetPosOffsetsParametric(point.X(), point.Y(), point.Z());
194  else
195  thePosOffsets.resize(3,0.0);
196  }
197 
198  return { thePosOffsets[0], thePosOffsets[1], thePosOffsets[2] };
199 }
200 
201 //----------------------------------------------------------------------------
203 std::vector<double> spacecharge::SpaceChargeStandard::GetPosOffsetsParametric(double xVal, double yVal, double zVal) const
204 {
205  std::vector<double> thePosOffsetsParametric;
206 
207  double xValNew = TransformX(xVal);
208  double yValNew = TransformY(yVal);
209  double zValNew = TransformZ(zVal);
210 
211  thePosOffsetsParametric.push_back(GetOnePosOffsetParametric(xValNew,yValNew,zValNew,"X"));
212  thePosOffsetsParametric.push_back(GetOnePosOffsetParametric(xValNew,yValNew,zValNew,"Y"));
213  thePosOffsetsParametric.push_back(GetOnePosOffsetParametric(xValNew,yValNew,zValNew,"Z"));
214 
215  return thePosOffsetsParametric;
216 }
217 
218 //----------------------------------------------------------------------------
221 double spacecharge::SpaceChargeStandard::GetOnePosOffsetParametric(double xValNew, double yValNew, double zValNew, std::string axis) const
222 {
223  double parA[6][7];
224  double parB[6];
225 
226  for(int j = 0; j < 6; j++)
227  {
228  for(int i = 0; i < 7; i++)
229  parA[j][i] = 0.0;
230 
231  parB[j] = 0.0;
232  }
233 
234  if(axis == "X")
235  {
236  for(int j = 0; j < 7; j++)
237  {
238  parA[0][j] = g1_x[j]->Eval(zValNew);
239  parA[1][j] = g2_x[j]->Eval(zValNew);
240  parA[2][j] = g3_x[j]->Eval(zValNew);
241  parA[3][j] = g4_x[j]->Eval(zValNew);
242  parA[4][j] = g5_x[j]->Eval(zValNew);
243  }
244 
245  f1_x->SetParameters(parA[0]);
246  f2_x->SetParameters(parA[1]);
247  f3_x->SetParameters(parA[2]);
248  f4_x->SetParameters(parA[3]);
249  f5_x->SetParameters(parA[4]);
250  }
251  else if(axis == "Y")
252  {
253  for(int j = 0; j < 6; j++)
254  {
255  parA[0][j] = g1_y[j]->Eval(zValNew);
256  parA[1][j] = g2_y[j]->Eval(zValNew);
257  parA[2][j] = g3_y[j]->Eval(zValNew);
258  parA[3][j] = g4_y[j]->Eval(zValNew);
259  parA[4][j] = g5_y[j]->Eval(zValNew);
260  parA[5][j] = g6_y[j]->Eval(zValNew);
261  }
262 
263  f1_y->SetParameters(parA[0]);
264  f2_y->SetParameters(parA[1]);
265  f3_y->SetParameters(parA[2]);
266  f4_y->SetParameters(parA[3]);
267  f5_y->SetParameters(parA[4]);
268  f6_y->SetParameters(parA[5]);
269  }
270  else if(axis == "Z")
271  {
272  for(int j = 0; j < 5; j++)
273  {
274  parA[0][j] = g1_z[j]->Eval(zValNew);
275  parA[1][j] = g2_z[j]->Eval(zValNew);
276  parA[2][j] = g3_z[j]->Eval(zValNew);
277  parA[3][j] = g4_z[j]->Eval(zValNew);
278  }
279 
280  f1_z->SetParameters(parA[0]);
281  f2_z->SetParameters(parA[1]);
282  f3_z->SetParameters(parA[2]);
283  f4_z->SetParameters(parA[3]);
284  }
285 
286  double aValNew;
287  double bValNew;
288 
289  if(axis == "Y")
290  {
291  aValNew = xValNew;
292  bValNew = yValNew;
293  }
294  else
295  {
296  aValNew = yValNew;
297  bValNew = xValNew;
298  }
299 
300  double offsetValNew = 0.0;
301  if(axis == "X")
302  {
303  parB[0] = f1_x->Eval(aValNew);
304  parB[1] = f2_x->Eval(aValNew);
305  parB[2] = f3_x->Eval(aValNew);
306  parB[3] = f4_x->Eval(aValNew);
307  parB[4] = f5_x->Eval(aValNew);
308 
309  fFinal_x->SetParameters(parB);
310  offsetValNew = 100.0*fFinal_x->Eval(bValNew);
311  }
312  else if(axis == "Y")
313  {
314  parB[0] = f1_y->Eval(aValNew);
315  parB[1] = f2_y->Eval(aValNew);
316  parB[2] = f3_y->Eval(aValNew);
317  parB[3] = f4_y->Eval(aValNew);
318  parB[4] = f5_y->Eval(aValNew);
319  parB[5] = f6_y->Eval(aValNew);
320 
321  fFinal_y->SetParameters(parB);
322  offsetValNew = 100.0*fFinal_y->Eval(bValNew);
323  }
324  else if(axis == "Z")
325  {
326  parB[0] = f1_z->Eval(aValNew);
327  parB[1] = f2_z->Eval(aValNew);
328  parB[2] = f3_z->Eval(aValNew);
329  parB[3] = f4_z->Eval(aValNew);
330 
331  fFinal_z->SetParameters(parB);
332  offsetValNew = 100.0*fFinal_z->Eval(bValNew);
333  }
334 
335  return offsetValNew;
336 }
337 
338 //----------------------------------------------------------------------------
342 {
343  std::vector<double> theEfieldOffsets;
344 
345  if(fRepresentationType == "Parametric")
346  theEfieldOffsets = GetEfieldOffsetsParametric(point.X(), point.Y(), point.Z());
347  else
348  theEfieldOffsets.resize(3,0.0);
349 
350  return { -theEfieldOffsets[0], -theEfieldOffsets[1], -theEfieldOffsets[2] };
351 }
352 
353 //----------------------------------------------------------------------------
355 std::vector<double> spacecharge::SpaceChargeStandard::GetEfieldOffsetsParametric(double xVal, double yVal, double zVal) const
356 {
357  std::vector<double> theEfieldOffsetsParametric;
358 
359  double xValNew = TransformX(xVal);
360  double yValNew = TransformY(yVal);
361  double zValNew = TransformZ(zVal);
362 
363  theEfieldOffsetsParametric.push_back(GetOneEfieldOffsetParametric(xValNew,yValNew,zValNew,"X"));
364  theEfieldOffsetsParametric.push_back(GetOneEfieldOffsetParametric(xValNew,yValNew,zValNew,"Y"));
365  theEfieldOffsetsParametric.push_back(GetOneEfieldOffsetParametric(xValNew,yValNew,zValNew,"Z"));
366 
367  return theEfieldOffsetsParametric;
368 }
369 
370 //----------------------------------------------------------------------------
373 double spacecharge::SpaceChargeStandard::GetOneEfieldOffsetParametric(double xValNew, double yValNew, double zValNew, std::string axis) const
374 {
375  double parA[6][7];
376  double parB[6];
377 
378  for(int j = 0; j < 6; j++)
379  {
380  for(int i = 0; i < 7; i++)
381  parA[j][i] = 0.0;
382 
383  parB[j] = 0.0;
384  }
385 
386  if(axis == "X")
387  {
388  for(int j = 0; j < 7; j++)
389  {
390  parA[0][j] = g1_Ex[j]->Eval(zValNew);
391  parA[1][j] = g2_Ex[j]->Eval(zValNew);
392  parA[2][j] = g3_Ex[j]->Eval(zValNew);
393  parA[3][j] = g4_Ex[j]->Eval(zValNew);
394  parA[4][j] = g5_Ex[j]->Eval(zValNew);
395  }
396 
397  f1_Ex->SetParameters(parA[0]);
398  f2_Ex->SetParameters(parA[1]);
399  f3_Ex->SetParameters(parA[2]);
400  f4_Ex->SetParameters(parA[3]);
401  f5_Ex->SetParameters(parA[4]);
402  }
403  else if(axis == "Y")
404  {
405  for(int j = 0; j < 6; j++)
406  {
407  parA[0][j] = g1_Ey[j]->Eval(zValNew);
408  parA[1][j] = g2_Ey[j]->Eval(zValNew);
409  parA[2][j] = g3_Ey[j]->Eval(zValNew);
410  parA[3][j] = g4_Ey[j]->Eval(zValNew);
411  parA[4][j] = g5_Ey[j]->Eval(zValNew);
412  parA[5][j] = g6_Ey[j]->Eval(zValNew);
413  }
414 
415  f1_Ey->SetParameters(parA[0]);
416  f2_Ey->SetParameters(parA[1]);
417  f3_Ey->SetParameters(parA[2]);
418  f4_Ey->SetParameters(parA[3]);
419  f5_Ey->SetParameters(parA[4]);
420  f6_Ey->SetParameters(parA[5]);
421  }
422  else if(axis == "Z")
423  {
424  for(int j = 0; j < 5; j++)
425  {
426  parA[0][j] = g1_Ez[j]->Eval(zValNew);
427  parA[1][j] = g2_Ez[j]->Eval(zValNew);
428  parA[2][j] = g3_Ez[j]->Eval(zValNew);
429  parA[3][j] = g4_Ez[j]->Eval(zValNew);
430  }
431 
432  f1_Ez->SetParameters(parA[0]);
433  f2_Ez->SetParameters(parA[1]);
434  f3_Ez->SetParameters(parA[2]);
435  f4_Ez->SetParameters(parA[3]);
436  }
437 
438  double aValNew;
439  double bValNew;
440 
441  if(axis == "Y")
442  {
443  aValNew = xValNew;
444  bValNew = yValNew;
445  }
446  else
447  {
448  aValNew = yValNew;
449  bValNew = xValNew;
450  }
451 
452  double offsetValNew = 0.0;
453  if(axis == "X")
454  {
455  parB[0] = f1_Ex->Eval(aValNew);
456  parB[1] = f2_Ex->Eval(aValNew);
457  parB[2] = f3_Ex->Eval(aValNew);
458  parB[3] = f4_Ex->Eval(aValNew);
459  parB[4] = f5_Ex->Eval(aValNew);
460 
461  fFinal_Ex->SetParameters(parB);
462  offsetValNew = fFinal_Ex->Eval(bValNew);
463  }
464  else if(axis == "Y")
465  {
466  parB[0] = f1_Ey->Eval(aValNew);
467  parB[1] = f2_Ey->Eval(aValNew);
468  parB[2] = f3_Ey->Eval(aValNew);
469  parB[3] = f4_Ey->Eval(aValNew);
470  parB[4] = f5_Ey->Eval(aValNew);
471  parB[5] = f6_Ey->Eval(aValNew);
472 
473  fFinal_Ey->SetParameters(parB);
474  offsetValNew = fFinal_Ey->Eval(bValNew);
475  }
476  else if(axis == "Z")
477  {
478  parB[0] = f1_Ez->Eval(aValNew);
479  parB[1] = f2_Ez->Eval(aValNew);
480  parB[2] = f3_Ez->Eval(aValNew);
481  parB[3] = f4_Ez->Eval(aValNew);
482 
483  fFinal_Ez->SetParameters(parB);
484  offsetValNew = fFinal_Ez->Eval(bValNew);
485  }
486 
487  return offsetValNew;
488 }
489 
490 //----------------------------------------------------------------------------
493 {
494  double xValNew;
495  xValNew = xVal/100.0;
496 
497  return xValNew;
498 }
499 
500 //----------------------------------------------------------------------------
503 {
504  double yValNew;
505  yValNew = yVal/100.0;
506 
507  return yValNew;
508 }
509 
510 //----------------------------------------------------------------------------
513 {
514  double zValNew;
515  zValNew = zVal/100.0;
516 
517  return zValNew;
518 }
519 
520 //----------------------------------------------------------------------------
522 bool spacecharge::SpaceChargeStandard::IsInsideBoundaries(double xVal, double yVal, double zVal) const
523 {
524  bool isInside = false;
525 
526  return isInside;
527 }
bool IsInsideBoundaries(double xVal, double yVal, double zVal) const
Check to see if point is inside boundaries of map - redefine this in experiment-specific implementati...
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:167
std::vector< double > GetEfieldOffsetsParametric(double xVal, double yVal, double zVal) const
Provides E field offsets using a parametric representation.
double GetOnePosOffsetParametric(double xVal, double yVal, double zVal, std::string axis) const
std::vector< double > GetPosOffsetsParametric(double xVal, double yVal, double zVal) const
Provides position offsets using a parametric representation.
double TransformY(double yVal) const
Transform Y to SCE Y coordinate - redefine this in experiment-specific implementation! ...
SpaceChargeStandard(fhicl::ParameterSet const &pset)
double TransformX(double xVal) const
Transform X to SCE X coordinate - redefine this in experiment-specific implementation! ...
geo::Vector_t GetEfieldOffsets(geo::Point_t const &point) const override
T get(std::string const &key) const
Definition: ParameterSet.h:231
double TransformZ(double zVal) const
Transform Z to SCE Z coordinate - redefine this in experiment-specific implementation! ...
bool has_key(std::string const &key) const
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
geo::Vector_t GetPosOffsets(geo::Point_t const &point) const override
bool EnableCorrSCE() const override
Return boolean indicating whether or not to apply SCE corrections.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:187
bool Configure(fhicl::ParameterSet const &pset)
double GetOneEfieldOffsetParametric(double xVal, double yVal, double zVal, std::string axis) const