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