EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PndTransMap.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PndTransMap.cxx
1 // -------------------------------------------------------------------------
2 // ----- PndTransMap source file -----
3 // ----- Created 29/01/07 by M. Al/Turany -----
4 // -------------------------------------------------------------------------
5 
6 #include "FairRun.h"
7 #include "FairRuntimeDb.h"
8 
9 
10 
11 #include "PndTransMap.h"
12 #include "PndTransPar.h"
13 #include "FairRunSim.h"
14 #include "TArrayF.h"
15 
16 
17 #include <iostream>
18 #include "stdlib.h"
19 
20 
21 
22 using namespace std;
23 // ------------- Default constructor ----------------------------------
24 PndTransMap::PndTransMap()
25  :PndFieldMap(), fHemiX(0), fHemiY(0), fBeamMom(0.)
26 {
27  fType = 4;
28 }
29 // ------------------------------------------------------------------------
30 
31 
32 
33 // ------------- Standard constructor ---------------------------------
34 PndTransMap::PndTransMap(const char* mapName,
35  const char* fileType, Double_t BeamMom)
36  : PndFieldMap(mapName, fileType), fHemiX(0), fHemiY(0),fBeamMom(BeamMom)
37 {
38  fType = 4;
39  TString Suffix="";
40  FairRunSim *fRun= FairRunSim::Instance();
41  if(fRun) fBeamMom= fRun->GetBeamMom();
42 
43  if(fBeamMom< 3)Suffix=".0150" ;
44  else if (fBeamMom< 6.0 && fBeamMom >= 3.0)Suffix=".0406";
45  else if (fBeamMom< 10.0 && fBeamMom >= 6.0 )Suffix=".0890" ;
46  else if (fBeamMom< 13.0 && fBeamMom >= 10.0)Suffix=".1191";
47  else if (fBeamMom> 13.0) Suffix=".1500";
48 
49 
50  TString NewName=mapName;
51  NewName=mapName+Suffix;
52  SetName(NewName.Data());
53  TString dir = getenv("VMCWORKDIR");
54  fFileName = dir + "/input/" + NewName;
55  if ( fileType[0] == 'R' ) fFileName += ".root";
56  else fFileName += ".dat";
57 
58 
59 
60 }
61 // ------------------------------------------------------------------------
62 
63 
64 
65 // ------------ Constructor from PndFieldPar --------------------------
66 PndTransMap::PndTransMap(PndTransPar* fieldPar)
67  : PndFieldMap(), fHemiX(0), fHemiY(0)
68 {
69  fType = 4;
70  fPosX = fPosY = fPosZ = 0.;
71  fXmin = fYmin = fZmin = 0.;
72  fXmax = fYmax = fZmax = 0.;
73  fXstep = fYstep = fZstep = 0.;
74  fNx = fNy = fNz = 0;
75  fScale = 1.;
76  funit = 10.0;
77  fBx = fBy = fBz = NULL;
78  if ( ! fieldPar ) {
79  cerr << "-W- PndTransMap::PndTransMap: empty parameter container!"
80  << endl;
81  fName = "";
82  fFileName = "";
83  fType = 1;
84  }
85  else {
86  fieldPar->MapName(fName);
87  fPosX = fieldPar->GetPositionX();
88  fPosY = fieldPar->GetPositionY();
89  fPosZ = fieldPar->GetPositionZ();
90  fScale = fieldPar->GetScale();
91  TString dir = getenv("VMCWORKDIR");
92  fFileName = dir + "/input/" + fName + ".root";
93  fType = fieldPar->GetType();
94  }
95 }
96 // ------------------------------------------------------------------------
97 
98 
99 
100 // ------------ Destructor --------------------------------------------
101 PndTransMap::~PndTransMap()
102 {
103  //printf("PndTransMap::~PndTransMap() \n ");
104 }
105 // ------------------------------------------------------------------------
106 void PndTransMap::GetBxyz(const Double_t point[3], Double_t* bField)
107 {
108  Double_t x =point[0];
109  Double_t y =point[1];
110  Double_t z =point[2];
111  Int_t ix = 0;
112  Int_t iy = 0;
113  Int_t iz = 0;
114  Double_t dx = 0.;
115  Double_t dy = 0.;
116  Double_t dz = 0.;
117 
118  if ( IsInside(x, y, z, ix, iy, iz, dx, dy, dz) ){
119  // Get Bx field values at grid cell corners
120  fHa[0][0][0] = fBx->At(ix *fNy*fNz + iy *fNz + iz);
121  fHa[1][0][0] = fBx->At((ix+1)*fNy*fNz + iy *fNz + iz);
122  fHa[0][1][0] = fBx->At(ix *fNy*fNz + (iy+1)*fNz + iz);
123  fHa[1][1][0] = fBx->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
124  fHa[0][0][1] = fBx->At(ix *fNy*fNz + iy *fNz + (iz+1));
125  fHa[1][0][1] = fBx->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
126  fHa[0][1][1] = fBx->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
127  fHa[1][1][1] = fBx->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
128 
129  //Bx is antisymtric in X
130  bField[0] =Interpolate(dx, dy, dz) * fHemiX ;
131 
132  // Get By field values at grid cell corners
133  fHa[0][0][0] = fBy->At(ix *fNy*fNz + iy *fNz + iz);
134  fHa[1][0][0] = fBy->At((ix+1)*fNy*fNz + iy *fNz + iz);
135  fHa[0][1][0] = fBy->At(ix *fNy*fNz + (iy+1)*fNz + iz);
136  fHa[1][1][0] = fBy->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
137  fHa[0][0][1] = fBy->At(ix *fNy*fNz + iy *fNz + (iz+1));
138  fHa[1][0][1] = fBy->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
139  fHa[0][1][1] = fBy->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
140  fHa[1][1][1] = fBy->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
141 
142  //By is symtric in X
143  bField[1] =Interpolate(dx, dy, dz);
144 
145  // Get Bz field values at grid cell corners
146  fHa[0][0][0] = fBz->At(ix *fNy*fNz + iy *fNz + iz);
147  fHa[1][0][0] = fBz->At((ix+1)*fNy*fNz + iy *fNz + iz);
148  fHa[0][1][0] = fBz->At(ix *fNy*fNz + (iy+1)*fNz + iz);
149  fHa[1][1][0] = fBz->At((ix+1)*fNy*fNz + (iy+1)*fNz + iz);
150  fHa[0][0][1] = fBz->At(ix *fNy*fNz + iy *fNz + (iz+1));
151  fHa[1][0][1] = fBz->At((ix+1)*fNy*fNz + iy *fNz + (iz+1));
152  fHa[0][1][1] = fBz->At(ix *fNy*fNz + (iy+1)*fNz + (iz+1));
153  fHa[1][1][1] = fBz->At((ix+1)*fNy*fNz + (iy+1)*fNz + (iz+1));
154 
155  // Return interpolated field value
156  //Bz is symtric in X
157  bField[2]=Interpolate(dx, dy, dz);
158 
159  }else{
160  bField[0]=0;
161  bField[1]=0;
162  bField[2]=0;
163  }
164 
165 }
166 
167 
168 
169 // ----------- Check whether a point is inside the map ----------------
170 Bool_t PndTransMap::IsInside(Double_t x, Double_t y, Double_t z,
171  Int_t& ix, Int_t& iy, Int_t& iz,
172  Double_t& dx, Double_t& dy,
173  Double_t& dz) {
174 
175  // --- Transform into local coordinate system
176  Double_t xl = x - fPosX;
177  Double_t yl = y - fPosY;
178  Double_t zl = z - fPosZ;
179 
180  // --- Reflect w.r.t. symmetry axes
181  fHemiX = fHemiY = 1.;
182  if ( xl < 0. ) {
183  fHemiX = -1.;
184  xl = -1. * xl;
185  }
186 
187  // --- Check for being outside the map range
188  if ( ! ( xl >= fXmin && xl <= fXmax && yl >= fYmin && yl <= fYmax &&
189  zl >= fZmin && zl <= fZmax ) ) {
190  ix = iy = iz = 0;
191  dx = dy = dz = 0.;
192  return kFALSE;
193  }
194 
195  // --- Determine grid cell
196  ix = Int_t( (xl-fXmin) / fXstep );
197  iy = Int_t( (yl-fYmin) / fYstep );
198  iz = Int_t( (zl-fZmin) / fZstep );
199 
200 
201  // Relative distance from grid point (in units of cell size)
202  dx = (xl-fXmin) / fXstep - Double_t(ix);
203  dy = (yl-fYmin) / fYstep - Double_t(iy);
204  dz = (zl-fZmin) / fZstep - Double_t(iz);
205 
206  return kTRUE;
207 
208 }
209 // ------------------------------------------------------------------------
210 void PndTransMap::FillParContainer()
211 {
212  TString MapName=GetName();
213  // cout << "PndConstField::FillParContainer() " << endl;
214  FairRun *fRun=FairRun::Instance();
215  FairRuntimeDb *rtdb=fRun->GetRuntimeDb();
216  Bool_t kParameterMerged=kTRUE;
217  PndTransPar* Par = (PndTransPar*) rtdb->getContainer("PndTransPar");
218  Par->SetParameters(this);
219  Par->setInputVersion(fRun->GetRunId(),1);
220  Par->setChanged();
221 
222 }
223 
224 
225 
226 
227 ClassImp(PndTransMap)