ATLAS Offline Software
Loading...
Searching...
No Matches
ForwardRegionFieldSvc.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "ForwardRegionFieldSvc.h"
6
7#include <fstream>
8
9#include "CLHEP/Geometry/Transform3D.h"
10
11#include "GaudiKernel/IIncidentSvc.h"
12
13#include "PathResolver/PathResolver.h"
14
15//comments kept as a reference
16// static ForwardRegionMgField q1("Q1",ForwardRegionField::Q1);
17// static ForwardRegionMgField q2("Q2",ForwardRegionField::Q2);
18// static ForwardRegionMgField q3("Q3",ForwardRegionField::Q3);
19// static ForwardRegionMgField d1("D1",ForwardRegionField::D1);
20// static ForwardRegionMgField d2("D2",ForwardRegionField::D2);
21// static ForwardRegionMgField q4("Q4",ForwardRegionField::Q4);
22// static ForwardRegionMgField q5("Q5",ForwardRegionField::Q5);
23// static ForwardRegionMgField q6("Q6",ForwardRegionField::Q6);
24// static ForwardRegionMgField q7("Q7",ForwardRegionField::Q7);
25// static ForwardRegionMgField q1hkick("Q1HKick",ForwardRegionField::Q1HKick);
26// static ForwardRegionMgField q1vkick("Q1VKick",ForwardRegionField::Q1VKick);
27// static ForwardRegionMgField q2hkick("Q2HKick",ForwardRegionField::Q2HKick);
28// static ForwardRegionMgField q2vkick("Q2VKick",ForwardRegionField::Q2VKick);
29// static ForwardRegionMgField q3hkick("Q3HKick",ForwardRegionField::Q3HKick);
30// static ForwardRegionMgField q3vkick("Q3VKick",ForwardRegionField::Q3VKick);
31// static ForwardRegionMgField q4vkickA("Q4VKickA",ForwardRegionField::Q4VKickA);
32// static ForwardRegionMgField q4hkick("Q4HKick",ForwardRegionField::Q4HKick);
33// static ForwardRegionMgField q4vkickB("Q4VKickB",ForwardRegionField::Q4VKickB);
34// static ForwardRegionMgField q5hkick("Q5HKick",ForwardRegionField::Q5HKick);
35// static ForwardRegionMgField q6vkick("Q6VKick",ForwardRegionField::Q6VKick);
36
37MagField::ForwardRegionFieldSvc::ForwardRegionFieldSvc(const std::string& name,ISvcLocator* svc) :
38 base_class(name,svc),
39 m_magnet(-1),
40 m_magDataType(0),
41 m_MQXA_DataFile(""), //"MQXA_NOMINAL.dat" if name = Q1 or Q3
42 m_MQXB_DataFile(""), //"MQXB_NOMINAL.dat" if name = Q2
43 m_properties("ForwardRegionProperties")
44{
45
46 declareProperty("Magnet", m_magnet, "integer property");
47 declareProperty("MQXA_DataFile", m_MQXA_DataFile, "");
48 declareProperty("MQXB_DataFile", m_MQXB_DataFile, "");
49 declareProperty("ForwardRegionProperties", m_properties);
50 m_Config.clear();
51
52 if(!m_MQXA_DataFile.empty())
53 {
54 const int fileMQXACols(7); //TODO make this a property
55
56 // load field map MQXA
57 std::vector<std::vector<std::string> > loadedDataFile = loadDataFile(m_MQXA_DataFile.c_str(),fileMQXACols);
58
59 ATH_MSG_DEBUG("Field map MQXA loaded");
60
61 // initialize magnetic induction mesh
62 for(int j=0; j < s_colsMQXA; j++)
63 {
64 for(int i=0; i < s_rowsMQXA; i++)
65 {
66 m_magIndMQXA[i][j][0] = atof(loadedDataFile[i+j*s_colsMQXA][4].c_str());
67 m_magIndMQXA[i][j][1] = atof(loadedDataFile[i+j*s_colsMQXA][5].c_str());
68 }
69 }
70 }
71
72 if(!m_MQXB_DataFile.empty())
73 {
74 const int fileMQXBCols(7); //TODO make this a property
75
76 // load field map MQXB
77 std::vector<std::vector<std::string> > loadedDataFileB = loadDataFile(m_MQXB_DataFile.c_str(),fileMQXBCols);
78
79 ATH_MSG_DEBUG("Field map MQXB loaded");
80
81 // initialize magnetic induction mesh
82 for(int j=0; j < s_colsMQXB; j++)
83 {
84 for(int i=0; i < s_rowsMQXB; i++)
85 {
86 m_magIndMQXB[i][j][0] = atof(loadedDataFileB[i+j*s_colsMQXB][4].c_str());
87 m_magIndMQXB[i][j][1] = atof(loadedDataFileB[i+j*s_colsMQXB][5].c_str());
88 }
89 }
90 }
91
92
93}
94
95StatusCode MagField::ForwardRegionFieldSvc::initialize()
96{
97 ATH_CHECK(m_properties.retrieve());
98
99 ATH_MSG_INFO("Postponing magnet strengths initialization of " << this->name() << " till the begin of run");
100 ServiceHandle<IIncidentSvc> incidentSvc("IncidentSvc", this->name());
101 ATH_CHECK(incidentSvc.retrieve());
102 incidentSvc->addListener( this, IncidentType::BeginRun );
103 ATH_MSG_INFO("Added listener to BeginRun incident");
104 return StatusCode::SUCCESS;
105}
106
107// Handle incident function - if BeginRun happens, initialize mag. fields
108void MagField::ForwardRegionFieldSvc::handle(const Incident& runIncident)
109{
110 ATH_MSG_INFO("handling incidents ..."); //FIXME drop to DEBUG level
111 if (runIncident.type() == IncidentType::BeginRun)
112 {
113 InitMagData();
114 }
115 ATH_MSG_INFO("BeginRun incident handled");
116}
117
118
123void MagField::ForwardRegionFieldSvc::getField(const double *xyz, double *bxyz, double*) const
124{
125 G4ThreeVector f = this->FieldValue(G4ThreeVector(xyz[0],xyz[1],xyz[2]));
126 bxyz[0] = f[0];
127 bxyz[1] = f[1];
128 bxyz[2] = f[2];
129}
130
135void MagField::ForwardRegionFieldSvc::getFieldZR(const double*, double*, double*) const
136{
137 throw std::logic_error("MagField::ForwardRegionFieldSvc::getFieldZR should not be called"); //FIXME not supported yet.
138}
139
140G4ThreeVector MagField::ForwardRegionFieldSvc::FieldValue(G4ThreeVector Point) const
141{
142 int beam;
143 if(Point[2] < 0) beam = 1;
144 else beam = 2;
145
146
147 HepGeom::Point3D<double> pointMagStart;
148 HepGeom::Point3D<double> pointMagEnd;
149 double rotZ = 0;
150
151 getMagnetTransformParams(beam, m_magnet, Point, pointMagStart, pointMagEnd, rotZ);
152
153 G4ThreeVector field;
154
155 // are points defined (non-zero)? If not, only rotation around Z will be applied (if any)
156 bool pointsNotDefined = (pointMagStart.distance2() == 0 || pointMagEnd.distance2() == 0);
157
158 // calculate x and y shifts -- for quadrupoles only, no effect in dipoles (magnet borders to be solved in GeoModel)
159 double xShift = pointsNotDefined ? getMagXOff(m_magnet) : (pointMagStart[0]+pointMagEnd[0])*0.5;
160 double yShift = pointsNotDefined ? getMagYOff(m_magnet) : (pointMagStart[1]+pointMagEnd[1])*0.5;
161
162 if(m_magnet <= s_Q3){ // inner triplet
163 if(m_Config.bUseFLUKAMapsForInnerTriplet)
164 field = getMagInd(Point,m_magnet, beam);
165 else {
166 double gradB = getMag(m_magnet,beam)*CLHEP::tesla/CLHEP::m*pow(-1.0,beam);
167 field = G4ThreeVector(gradB*(Point[1]-yShift),gradB*(Point[0]-xShift),0);
168 }
169 }
170 else if(m_magnet <= s_D2) // dipoles
171 field = G4ThreeVector(0,getMag(m_magnet,beam)*CLHEP::tesla,0);
172 else if(m_magnet <= s_Q7) // other quadrupoles
173 {
174 double gradB = getMag(m_magnet,beam)*CLHEP::tesla/CLHEP::m*pow(-1.0,beam);
175 field = G4ThreeVector(gradB*(Point[1]-yShift),gradB*(Point[0]-xShift),0);
176 }
177 else // kickers
178 return getKick(beam);
179
180 // No points defined, check if there is rotation around magnet axis, else return unchanged field
181 if(pointsNotDefined){
182 if(rotZ == 0) return field;
183 return HepGeom::RotateZ3D(rotZ)*(HepGeom::Vector3D<double>)field;
184 }
185
186 // Determine rotation from start and end points, shift is already taken care of
187 HepGeom::Point3D<double> pointMagEndNoShift(pointMagEnd[0]-xShift, pointMagEnd[1]-yShift, pointMagEnd[2]);
188
189 HepGeom::Point3D<double> pointMagRotCenter(0,0,(pointMagStart[2]+pointMagEnd[2])/2);
190
191 HepGeom::Vector3D<double> vecNoRot(0,0,pointMagEnd[2]-pointMagRotCenter[2]);
192 HepGeom::Vector3D<double> vecRot = pointMagEndNoShift - pointMagRotCenter;
193
194 HepGeom::Vector3D<double> vecRotAxis = vecNoRot.cross(vecRot);
195 double angle = vecNoRot.angle(vecRot);
196
197 HepGeom::Transform3D rotateField = HepGeom::Rotate3D(angle, vecRotAxis);
198
199 field = rotateField*HepGeom::RotateZ3D(rotZ)*(HepGeom::Vector3D<double>)field;
200 return field;
201}
202
203
204
205
206
207
208
209
210
211
212// The main initialization - initialization of fields based on twiss files or magnets.dat
213void MagField::ForwardRegionFieldSvc::InitMagData()
214{
215 ATH_MSG_INFO("Initializing magnetic field");
216
217 m_Config = *(m_properties->getConf());
218
219 if(m_Config.twissFileB1 == "" || m_Config.twissFileB2 == "" || m_Config.momentum == 0){
220 m_magDataType = 0;
221
222 ATH_MSG_INFO("Using magnets.dat as the field settings source");
223
224 m_magnets = loadDataFile("ForwardRegionMgField/magnets.dat" ,5);
225 }
226 else
227 {
228 m_magDataType = 1;
229
230 ATH_MSG_INFO("Using twiss files (" << m_Config.twissFileB1 << ", " << m_Config.twissFileB2 << ") as the field settings source");
231
232 std::string headerB1;
233 std::string headerB2;
234 // Load twiss files and calculate field values for magnets
235 std::vector<std::vector<std::string> > loadedTwissFileB1 = loadDataFileNLines(m_Config.twissFileB1.c_str(),30,200,headerB1);
236 std::vector<std::vector<std::string> > loadedTwissFileB2 = loadDataFileNLines(m_Config.twissFileB2.c_str(),30,200,headerB2);
237
238 // Beam 1
239 InitMagDataFromTwiss(loadedTwissFileB1, 1, m_Config.momentum);
240
241 // Beam 2
242 InitMagDataFromTwiss(loadedTwissFileB2, 2, m_Config.momentum);
243
244 //writeOutTwiss(loadedTwissFileB1,1,headerB1);
245 //writeOutTwiss(loadedTwissFileB2,2,headerB2);
246 ATH_MSG_INFO("Field initialized.");
247 }
248}
249
250void MagField::ForwardRegionFieldSvc::InitMagDataFromTwiss(const std::vector<std::vector<std::string> >& loadedTwissFile, int beam, double momentum)
251{
252 int textIndex = 0;
253 int lengthIndex = 4;
254 int hkick = 5;
255 int vkick = 6;
256 int k0LIndex = 7;
257 int k1LIndex = 8;
258
259 bool dipole = false;
260 bool quadrupole = false;
261 bool kicker = false;
262 double length;
263
264 // init offsets (Q4--Q7 have dx=-97 mm, others 0), do not apply to kickers
265 if(m_magnet <= s_Q7) {
266 if(m_magnet >= s_Q4)
267 m_magData[2] = -97*CLHEP::mm;
268 else
269 m_magData[2] = 0;
270
271 m_magData[3] = 0;
272 }
273
274 for(int i=0; i < 150; i++)
275 {
276 // dipoles
277 dipole = (m_magnet == s_D1 && (loadedTwissFile[i][textIndex] == "\"MBXW.A4R1\"" || loadedTwissFile[i][textIndex] == "\"MBXW.A4L1\""))
278 || (m_magnet == s_D2 && (loadedTwissFile[i][textIndex] == "\"MBRC.4R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MBRC.4L1.B2\""));
279 if(dipole){
280 m_magData[beam-1] = kLToB(atof(loadedTwissFile[i][k0LIndex].c_str()), atof(loadedTwissFile[i][lengthIndex].c_str()), momentum);
281 return;
282 }
283
284 //quadrupoles
285 quadrupole = (m_magnet == s_Q1 && (loadedTwissFile[i][textIndex] == "\"MQXA.1R1\"" || loadedTwissFile[i][textIndex] == "\"MQXA.1L1\""))
286 || (m_magnet == s_Q2 && (loadedTwissFile[i][textIndex] == "\"MQXB.A2R1\"" || loadedTwissFile[i][textIndex] == "\"MQXB.A2L1\""))
287 || (m_magnet == s_Q3 && (loadedTwissFile[i][textIndex] == "\"MQXA.3R1\"" || loadedTwissFile[i][textIndex] == "\"MQXA.3L1\""))
288 || (m_magnet == s_Q4 && (loadedTwissFile[i][textIndex] == "\"MQY.4R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MQY.4L1.B2\""))
289 || (m_magnet == s_Q5 && (loadedTwissFile[i][textIndex] == "\"MQML.5R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MQML.5L1.B2\""))
290 || (m_magnet == s_Q6 && (loadedTwissFile[i][textIndex] == "\"MQML.6R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MQML.6L1.B2\""))
291 || (m_magnet == s_Q7 && (loadedTwissFile[i][textIndex] == "\"MQM.A7R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MQM.A7L1.B2\""));
292 if(quadrupole){
293 m_magData[beam-1] = kLToB(atof(loadedTwissFile[i][k1LIndex].c_str()), atof(loadedTwissFile[i][lengthIndex].c_str()), momentum);
294 return;
295 }
296
297 // kickers
298 kicker = (m_magnet == s_Q1HKick && (loadedTwissFile[i][textIndex] == "\"MCBXH.1R1\"" || loadedTwissFile[i][textIndex] == "\"MCBXH.1L1\""))
299 || (m_magnet == s_Q1VKick && (loadedTwissFile[i][textIndex] == "\"MCBXV.1R1\"" || loadedTwissFile[i][textIndex] == "\"MCBXV.1L1\""))
300 || (m_magnet == s_Q2HKick && (loadedTwissFile[i][textIndex] == "\"MCBXH.2R1\"" || loadedTwissFile[i][textIndex] == "\"MCBXH.2L1\""))
301 || (m_magnet == s_Q2VKick && (loadedTwissFile[i][textIndex] == "\"MCBXV.2R1\"" || loadedTwissFile[i][textIndex] == "\"MCBXV.2L1\""))
302 || (m_magnet == s_Q3HKick && (loadedTwissFile[i][textIndex] == "\"MCBXH.3R1\"" || loadedTwissFile[i][textIndex] == "\"MCBXH.3L1\""))
303 || (m_magnet == s_Q3VKick && (loadedTwissFile[i][textIndex] == "\"MCBXV.3R1\"" || loadedTwissFile[i][textIndex] == "\"MCBXV.3L1\""))
304 || (m_magnet == s_Q4VKickA && (loadedTwissFile[i][textIndex] == "\"MCBYV.A4R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MCBYV.A4L1.B2\""))
305 || (m_magnet == s_Q4HKick && (loadedTwissFile[i][textIndex] == "\"MCBYH.4R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MCBYH.4L1.B2\""))
306 || (m_magnet == s_Q4VKickB && (loadedTwissFile[i][textIndex] == "\"MCBYV.B4R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MCBYV.B4L1.B2\""))
307 || (m_magnet == s_Q5HKick && (loadedTwissFile[i][textIndex] == "\"MCBCH.5R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MCBCH.5L1.B2\""))
308 || (m_magnet == s_Q6VKick && (loadedTwissFile[i][textIndex] == "\"MCBCV.6R1.B1\"" || loadedTwissFile[i][textIndex] == "\"MCBCV.6L1.B2\""));
309
310 if(kicker){
311 // length of s_Q1--s_Q3 kickers in twiss files is zero -> it is set to 0.001 m (also the length of corresponding GeoModel volumes)
312 length = atof(loadedTwissFile[i][lengthIndex].c_str()) ? atof(loadedTwissFile[i][lengthIndex].c_str()) : 0.001;
313 // calculation of the B from deflection angle --- beam going down at the IP
314 m_magData[beam-1] = pow(-1.0, beam+1)*kLToB(atof(loadedTwissFile[i][vkick].c_str()), length, momentum);
315 m_magData[beam+1] = kLToB(atof(loadedTwissFile[i][hkick].c_str()), length, momentum);
316 return;
317 }
318 }
319}
320
321
323// get magnetic induction vector in certain point inside MQXA (MQXAB = 0) or MQXB (MQXAB = 1) (uses bilinear interpolation), q=0 - s_Q1, q=1 - Q2a, q=2 - Q2b, q=3 - s_Q3
324G4ThreeVector MagField::ForwardRegionFieldSvc::getMagInd(G4ThreeVector Point, int q, int beam) const
325{
326 int MQXAB = 0;
327 if(q == 1) MQXAB = 1;
328
329 //magnet gradient and offsets
330 double gradMQX = getMag(q,beam);
331 double xOffMQX = getMagXOff(q);
332 double yOffMQX = getMagYOff(q);
333
334 //mapfile parameters
335 double minMQXA = -24*CLHEP::cm;
336 double maxMQXA = 24*CLHEP::cm;
337 double minMQXB = -26*CLHEP::cm;
338 double maxMQXB = 26*CLHEP::cm;
339
340 double x,y,x1,y1,x2,y2,xstep,ystep;
341 int i1,i2,j1,j2;
342
343 G4ThreeVector B(0,0,0);
344 if(MQXAB == 0) // field in MQXA
345 {
346 x = Point[0] - xOffMQX;
347 y = Point[1] - yOffMQX;
348
349 xstep = (maxMQXA - minMQXA)/(s_rowsMQXA-1);
350 ystep = (maxMQXA - minMQXA)/(s_colsMQXA-1);
351
352
353 i2 = ceil(x/xstep)+ceil(s_rowsMQXA/2.);
354 j2 = ceil(y/ystep)+ceil(s_colsMQXA/2.);
355 i1 = i2 - 1;
356 j1 = j2 - 1;
357 x2 = (i2-ceil(s_rowsMQXA/2.))*xstep;
358 y2 = (j2-ceil(s_colsMQXA/2.))*ystep;
359 x1 = x2 - xstep;
360 y1 = y2 - ystep;
361
362 if(x < maxMQXA && x > minMQXA && y < maxMQXA && y > minMQXA) // field in MQXA
363 {
364 B[0] = (m_magIndMQXA[i1][j1][0]*(x2-x)*(y2-y) + m_magIndMQXA[i2][j1][0]*(x-x1)*(y2-y) + m_magIndMQXA[i1][j2][0]*(x2-x)*(y-y1) + m_magIndMQXA[i2][j2][0]*(x-x1)*(y-y1))/xstep/ystep*CLHEP::tesla;
365 B[1] = (m_magIndMQXA[i1][j1][1]*(x2-x)*(y2-y) + m_magIndMQXA[i2][j1][1]*(x-x1)*(y2-y) + m_magIndMQXA[i1][j2][1]*(x2-x)*(y-y1) + m_magIndMQXA[i2][j2][1]*(x-x1)*(y-y1))/xstep/ystep*CLHEP::tesla;
366 B[0] = B[0]*pow(-1.0,beam)*(gradMQX/224.29);
367 B[1] = B[1]*pow(-1.0,beam)*(gradMQX/224.29);
368 }
369 }
370 if(MQXAB == 1) // field in MQXB
371 {
372 x = Point[0] - xOffMQX;
373 y = Point[1] - yOffMQX;
374
375 xstep = (maxMQXB - minMQXB)/(s_rowsMQXB-1);
376 ystep = (maxMQXB - minMQXB)/(s_colsMQXB-1);
377
378 i2 = ceil(x/xstep)+ceil(s_rowsMQXB/2);
379 j2 = ceil(y/ystep)+ceil(s_colsMQXB/2);
380 i1 = i2 - 1;
381 j1 = j2 - 1;
382 x2 = (i2-ceil(s_rowsMQXB/2))*xstep;
383 y2 = (j2-ceil(s_colsMQXB/2))*ystep;
384 x1 = x2 - xstep;
385 y1 = y2 - ystep;
386
387 if(x < maxMQXB && x > minMQXB && y < maxMQXB && y > minMQXB) // field in MQXB
388 {
389 B[0] = (m_magIndMQXB[i1][j1][0]*(x2-x)*(y2-y) + m_magIndMQXB[i2][j1][0]*(x-x1)*(y2-y) + m_magIndMQXB[i1][j2][0]*(x2-x)*(y-y1) + m_magIndMQXB[i2][j2][0]*(x-x1)*(y-y1))/xstep/ystep*CLHEP::tesla;
390 B[1] = (m_magIndMQXB[i1][j1][1]*(x2-x)*(y2-y) + m_magIndMQXB[i2][j1][1]*(x-x1)*(y2-y) + m_magIndMQXB[i1][j2][1]*(x2-x)*(y-y1) + m_magIndMQXB[i2][j2][1]*(x-x1)*(y-y1))/xstep/ystep*CLHEP::tesla;
391 B[0] = B[0]*pow(-1.0,beam)*(gradMQX/216.1787104);
392 B[1] = B[1]*pow(-1.0,beam)*(gradMQX/216.1787104);
393 }
394 }
395 return B;
396}
397
398double MagField::ForwardRegionFieldSvc::getMag(int magnet, int beam) const
399{
400 if(m_magDataType == 1) return m_magData[beam-1];
401 return atof(m_magnets[magnet][beam].c_str());
402}
403
404G4ThreeVector MagField::ForwardRegionFieldSvc::getKick(int beam) const
405{
406 if(m_magDataType == 1) return G4ThreeVector(m_magData[beam-1]*CLHEP::tesla,m_magData[beam+1]*CLHEP::tesla,0);
407 return G4ThreeVector(0,0,0); // magnets.dat does not contain kickers
408}
409
410double MagField::ForwardRegionFieldSvc::getMagXOff(int magnet) const
411{
412 if(m_magDataType == 1) return m_magData[2];
413 return atof(m_magnets[magnet][3].c_str())*CLHEP::mm;
414}
415
416double MagField::ForwardRegionFieldSvc::getMagYOff(int magnet) const
417{
418 if(m_magDataType == 1) return m_magData[3];
419 return atof(m_magnets[magnet][4].c_str())*CLHEP::mm;
420}
421
422// Load data from file into 2D array of strings. Input is filename, wanted number of columns and max number of lines (in original file) to be loaded (N=0 means all)
423std::vector<std::vector<std::string> > MagField::ForwardRegionFieldSvc::loadDataFileNLines(const char * fileName, int cols, int N, std::string& header)
424{
425 std::vector<std::vector<std::string> > loadedData;
426
427 header = "";
428
429 MsgStream LogStream(Athena::getMessageSvc(), "MagField::ForwardRegionFieldSvc::loadDataFile()");
430
431 std::ifstream file (fileName);
432 if(!file){
433 std::string datapath = PathResolverFindDataFile(fileName);
434 LogStream << MSG::DEBUG << "File " << fileName << " not found in run directory, trying to load it from DATAPATH" << endmsg;
435 file.open(datapath.c_str());
436 }
437
438 if(!file)
439 LogStream << MSG::FATAL << "Unable to load " << fileName << endmsg;
440
441 if(file.is_open())
442 {
443 std::vector<std::string> row (cols);
444 char c;
445 char temp[1024];
446
447 bool unlimited = false;
448 if(N==0) unlimited = true;
449 int i = 0;
450 while(file.get(c) && (i < N || unlimited))
451 {
452 if(c != '@' && c != '#' && c != '*' && c != '$' && c != '%' && c!='\n')
453 {
454 file.unget();
455 for(int i = 0; i<cols; i++) // load desired columns
456 {
457 file >> row[i];
458 }
459 loadedData.push_back(row); //store them
460 file.ignore(1024,'\n'); // discard rest of line
461 }
462 else if(c == '@')
463 {
464 file.unget();
465 file.getline(temp,1024);
466 header.append(temp);
467 header.append("\n");
468 }
469 else if(c == '*')
470 {
471 file.unget();
472 file.getline(temp,1024);
473 header.append(temp);
474 header.append(" PosXStart");
475 header.append(" PosYStart");
476 header.append(" PosZStart");
477 header.append(" PosXEnd ");
478 header.append(" PosYEnd ");
479 header.append(" PosZEnd ");
480 header.append(" RotZ");
481 header.append("\n");
482 }
483 else if(c == '$')
484 {
485 file.unget();
486 file.getline(temp,1024);
487 header.append(temp);
488 header.append(" %le ");
489 header.append(" %le ");
490 header.append(" %le ");
491 header.append(" %le ");
492 header.append(" %le ");
493 header.append(" %le ");
494 header.append(" %le ");
495 header.append("\n");
496 }
497 else
498 file.ignore(1024, '\n'); // discard commented lines
499 i++;
500 }
501 LogStream << MSG::INFO << "File " << fileName << " succesfully loaded." << endmsg;
502 file.close();
503 }
504 return loadedData;
505}
506
507std::vector<std::vector<std::string> > MagField::ForwardRegionFieldSvc::loadDataFile(const char * fileName, int cols)
508{
509 std::string header;
510 return loadDataFileNLines(fileName, cols, 0, header);
511}
512
513// Calculate magnetic induction value / gradient from K*L values, same relation holds for induction value from deviation angle (hkick, vkick)
514double MagField::ForwardRegionFieldSvc::kLToB(double kL, double length, double momentum)
515{
516 return kL*momentum/length/0.299792458;
517}
518
519void MagField::ForwardRegionFieldSvc::writeOutTwiss(const std::vector<std::vector<std::string>>& loadedTwissFile, int beam, const std::string& header)
520{
521 MsgStream LogStream(Athena::getMessageSvc(), "MagField::ForwardRegionFieldSvc::writeOutTwiss()");
522
523
524 std::ostringstream fileName;
525 fileName << "alfaTwiss" << beam << ".txt";
526 std::ofstream file (fileName.str().c_str());
527 int magID;
528 HepGeom::Point3D<double> pointMagStart;
529 HepGeom::Point3D<double> pointMagEnd;
530 double rotZ = 0;
531 double PointZ;
532 G4ThreeVector Point;
533
534 if(!file)
535 LogStream << MSG::ERROR << "Unable to write to " << fileName.str() << endmsg;
536
537 if(file.is_open())
538 {
539 file << header;
540 int fileSize = loadedTwissFile.size();
541 int rowSize = 0;
542 for(int i=0; i < fileSize;i++)
543 {
544 rowSize = loadedTwissFile[i].size();
545 for(int j=0; j<rowSize; j++)
546 {
547 if(j < 3) file << std::left;
548 else file << std::right;
549 file << std::setw(20) << loadedTwissFile[i][j];
550 }
551 magID = getMagNumFromName(loadedTwissFile[i][0]);
552 if(magID >= 0){
553 PointZ = atof(loadedTwissFile[i][3].c_str())*CLHEP::m;
554 Point.set(0,0,PointZ);
555 getMagnetTransformParams(beam, magID, Point, pointMagStart, pointMagEnd, rotZ);
556 file << std::setw(20) << pointMagStart[0]/CLHEP::m;
557 file << std::setw(20) << pointMagStart[1]/CLHEP::m;
558 file << std::setw(20) << pointMagStart[2]/CLHEP::m;
559 file << std::setw(20) << pointMagEnd[0]/CLHEP::m;
560 file << std::setw(20) << pointMagEnd[1]/CLHEP::m;
561 file << std::setw(20) << pointMagEnd[2]/CLHEP::m;
562 file << std::setw(20) << rotZ;
563 }
564 else {
565 file << std::setw(20) << 0;
566 file << std::setw(20) << 0;
567 file << std::setw(20) << 0;
568 file << std::setw(20) << 0;
569 file << std::setw(20) << 0;
570 file << std::setw(20) << 0;
571 file << std::setw(20) << 0;
572 }
573
574 file << std::endl;
575 }
576 }
577 file.close();
578}
579
580int MagField::ForwardRegionFieldSvc::getMagNumFromName(const std::string& name) const
581{
582 // dipoles
583 if (name == "\"MBXW.A4R1\"" || name == "\"MBXW.A4L1\"") return s_D1;
584 if (name == "\"MBXW.B4R1\"" || name == "\"MBXW.B4L1\"") return s_D1;
585 if (name == "\"MBXW.C4R1\"" || name == "\"MBXW.C4L1\"") return s_D1;
586 if (name == "\"MBXW.D4R1\"" || name == "\"MBXW.D4L1\"") return s_D1;
587 if (name == "\"MBXW.E4R1\"" || name == "\"MBXW.E4L1\"") return s_D1;
588 if (name == "\"MBXW.F4R1\"" || name == "\"MBXW.F4L1\"") return s_D1;
589 if (name == "\"MBRC.4R1.B1\"" || name == "\"MBRC.4L1.B2\"") return s_D2;
590
591 //quadrupoles
592 if(name == "\"MQXA.1R1\"" || name == "\"MQXA.1L1\"") return s_Q1;
593 if(name == "\"MQXB.A2R1\"" || name == "\"MQXB.A2L1\"") return s_Q2;
594 if(name == "\"MQXB.B2R1\"" || name == "\"MQXB.B2L1\"") return s_Q2;
595 if(name == "\"MQXA.3R1\"" || name == "\"MQXA.3L1\"") return s_Q3;
596 if(name == "\"MQY.4R1.B1\"" || name == "\"MQY.4L1.B2\"") return s_Q4;
597 if(name == "\"MQML.5R1.B1\"" || name == "\"MQML.5L1.B2\"") return s_Q5;
598 if(name == "\"MQML.6R1.B1\"" || name == "\"MQML.6L1.B2\"") return s_Q6;
599 if(name == "\"MQM.A7R1.B1\"" || name == "\"MQM.A7L1.B2\"") return s_Q7;
600 if(name == "\"MQM.B7R1.B1\"" || name == "\"MQM.B7L1.B2\"") return s_Q7;
601
602 // kickers
603 if(name == "\"MCBXH.1R1\"" || name == "\"MCBXH.1L1\"") return s_Q1HKick;
604 if(name == "\"MCBXV.1R1\"" || name == "\"MCBXV.1L1\"") return s_Q1VKick;
605 if(name == "\"MCBXH.2R1\"" || name == "\"MCBXH.2L1\"") return s_Q2HKick;
606 if(name == "\"MCBXV.2R1\"" || name == "\"MCBXV.2L1\"") return s_Q2VKick;
607 if(name == "\"MCBXH.3R1\"" || name == "\"MCBXH.3L1\"") return s_Q3HKick;
608 if(name == "\"MCBXV.3R1\"" || name == "\"MCBXV.3L1\"") return s_Q3VKick;
609 if(name == "\"MCBYV.A4R1.B1\"" || name == "\"MCBYV.A4L1.B2\"") return s_Q4VKickA;
610 if(name == "\"MCBYH.4R1.B1\"" || name == "\"MCBYH.4L1.B2\"") return s_Q4HKick;
611 if(name == "\"MCBYV.B4R1.B1\"" || name == "\"MCBYV.B4L1.B2\"") return s_Q4VKickB;
612 if(name == "\"MCBCH.5R1.B1\"" || name == "\"MCBCH.5L1.B2\"") return s_Q5HKick;
613 if(name == "\"MCBCV.6R1.B1\"" || name == "\"MCBCV.6L1.B2\"") return s_Q6VKick;
614 return -1;
615}
616
617void MagField::ForwardRegionFieldSvc::getMagnetTransformParams(int beam, int magnet, G4ThreeVector Point, HepGeom::Point3D<double> &pointMagStart,HepGeom::Point3D<double> &pointMagEnd, double &rotZ) const
618{
619 // find out which magnet we are in and get corresponding displacements
620 switch(magnet){
621 case s_Q1:
622 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ1Start[(beam-1)*3],m_Config.pointQ1Start[(beam-1)*3+1],m_Config.pointQ1Start[(beam-1)*3+2]);
623 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ1End[(beam-1)*3],m_Config.pointQ1End[(beam-1)*3+1],m_Config.pointQ1End[(beam-1)*3+2]);
624 rotZ = m_Config.fQ1RotZ[beam-1];
625 break;
626 case s_Q2:
627 if(std::abs(Point[2]) < 38*CLHEP::m){
628 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ2aStart[(beam-1)*3],m_Config.pointQ2aStart[(beam-1)*3+1],m_Config.pointQ2aStart[(beam-1)*3+2]);
629 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ2aEnd[(beam-1)*3],m_Config.pointQ2aEnd[(beam-1)*3+1],m_Config.pointQ2aEnd[(beam-1)*3+2]);
630 rotZ = m_Config.fQ2aRotZ[beam-1];
631 }
632 else
633 {
634 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ2bStart[(beam-1)*3],m_Config.pointQ2bStart[(beam-1)*3+1],m_Config.pointQ2bStart[(beam-1)*3+2]);
635 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ2bEnd[(beam-1)*3],m_Config.pointQ2bEnd[(beam-1)*3+1],m_Config.pointQ2bEnd[(beam-1)*3+2]);
636 rotZ = m_Config.fQ2bRotZ[beam-1];
637 }
638 break;
639 case s_Q3:
640 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ3Start[(beam-1)*3],m_Config.pointQ3Start[(beam-1)*3+1],m_Config.pointQ3Start[(beam-1)*3+2]);
641 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ3End[(beam-1)*3],m_Config.pointQ3End[(beam-1)*3+1],m_Config.pointQ3End[(beam-1)*3+2]);
642 rotZ = m_Config.fQ3RotZ[beam-1];
643 break;
644 case s_Q4:
645 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ4Start[(beam-1)*3],m_Config.pointQ4Start[(beam-1)*3+1],m_Config.pointQ4Start[(beam-1)*3+2]);
646 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ4End[(beam-1)*3],m_Config.pointQ4End[(beam-1)*3+1],m_Config.pointQ4End[(beam-1)*3+2]);
647 rotZ = m_Config.fQ4RotZ[beam-1];
648 break;
649 case s_Q5:
650 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ5Start[(beam-1)*3],m_Config.pointQ5Start[(beam-1)*3+1],m_Config.pointQ5Start[(beam-1)*3+2]);
651 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ5End[(beam-1)*3],m_Config.pointQ5End[(beam-1)*3+1],m_Config.pointQ5End[(beam-1)*3+2]);
652 rotZ = m_Config.fQ5RotZ[beam-1];
653 break;
654 case s_Q6:
655 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ6Start[(beam-1)*3],m_Config.pointQ6Start[(beam-1)*3+1],m_Config.pointQ6Start[(beam-1)*3+2]);
656 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ6End[(beam-1)*3],m_Config.pointQ6End[(beam-1)*3+1],m_Config.pointQ6End[(beam-1)*3+2]);
657 rotZ = m_Config.fQ6RotZ[beam-1];
658 break;
659 case s_Q7:
660 if(std::abs(Point[2]) < 263.5*CLHEP::m){
661 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ7aStart[(beam-1)*3],m_Config.pointQ7aStart[(beam-1)*3+1],m_Config.pointQ7aStart[(beam-1)*3+2]);
662 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ7aEnd[(beam-1)*3],m_Config.pointQ7aEnd[(beam-1)*3+1],m_Config.pointQ7aEnd[(beam-1)*3+2]);
663 rotZ = m_Config.fQ7aRotZ[beam-1];
664 }
665 else
666 {
667 pointMagStart = HepGeom::Point3D<double>(m_Config.pointQ7bStart[(beam-1)*3],m_Config.pointQ7bStart[(beam-1)*3+1],m_Config.pointQ7bStart[(beam-1)*3+2]);
668 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointQ7bEnd[(beam-1)*3],m_Config.pointQ7bEnd[(beam-1)*3+1],m_Config.pointQ7bEnd[(beam-1)*3+2]);
669 rotZ = m_Config.fQ7bRotZ[beam-1];
670 }
671 break;
672 case s_D1:
673 if(std::abs(Point[2]) < 63.5*CLHEP::m){
674 pointMagStart = HepGeom::Point3D<double>(m_Config.pointD1aStart[(beam-1)*3],m_Config.pointD1aStart[(beam-1)*3+1],m_Config.pointD1aStart[(beam-1)*3+2]);
675 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointD1aEnd[(beam-1)*3],m_Config.pointD1aEnd[(beam-1)*3+1],m_Config.pointD1aEnd[(beam-1)*3+2]);
676 rotZ = m_Config.fD1aRotZ[beam-1];
677 }
678 else if(std::abs(Point[2]) < 67.5*CLHEP::m){
679 pointMagStart = HepGeom::Point3D<double>(m_Config.pointD1bStart[(beam-1)*3],m_Config.pointD1bStart[(beam-1)*3+1],m_Config.pointD1bStart[(beam-1)*3+2]);
680 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointD1bEnd[(beam-1)*3],m_Config.pointD1bEnd[(beam-1)*3+1],m_Config.pointD1bEnd[(beam-1)*3+2]);
681 rotZ = m_Config.fD1bRotZ[beam-1];
682 }
683
684 else if(std::abs(Point[2]) < 72*CLHEP::m){
685 pointMagStart = HepGeom::Point3D<double>(m_Config.pointD1cStart[(beam-1)*3],m_Config.pointD1cStart[(beam-1)*3+1],m_Config.pointD1cStart[(beam-1)*3+2]);
686 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointD1cEnd[(beam-1)*3],m_Config.pointD1cEnd[(beam-1)*3+1],m_Config.pointD1cEnd[(beam-1)*3+2]);
687 rotZ = m_Config.fD1cRotZ[beam-1];
688 }
689
690 else if(std::abs(Point[2]) < 76*CLHEP::m){
691 pointMagStart = HepGeom::Point3D<double>(m_Config.pointD1dStart[(beam-1)*3],m_Config.pointD1dStart[(beam-1)*3+1],m_Config.pointD1dStart[(beam-1)*3+2]);
692 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointD1dEnd[(beam-1)*3],m_Config.pointD1dEnd[(beam-1)*3+1],m_Config.pointD1dEnd[(beam-1)*3+2]);
693 rotZ = m_Config.fD1dRotZ[beam-1];
694 }
695
696 else if(std::abs(Point[2]) < 80.5*CLHEP::m){
697 pointMagStart = HepGeom::Point3D<double>(m_Config.pointD1eStart[(beam-1)*3],m_Config.pointD1eStart[(beam-1)*3+1],m_Config.pointD1eStart[(beam-1)*3+2]);
698 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointD1eEnd[(beam-1)*3],m_Config.pointD1eEnd[(beam-1)*3+1],m_Config.pointD1eEnd[(beam-1)*3+2]);
699 rotZ = m_Config.fD1eRotZ[beam-1];
700 }
701
702 else {
703 pointMagStart = HepGeom::Point3D<double>(m_Config.pointD1fStart[(beam-1)*3],m_Config.pointD1fStart[(beam-1)*3+1],m_Config.pointD1fStart[(beam-1)*3+2]);
704 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointD1fEnd[(beam-1)*3],m_Config.pointD1fEnd[(beam-1)*3+1],m_Config.pointD1fEnd[(beam-1)*3+2]);
705 rotZ = m_Config.fD1fRotZ[beam-1];
706 }
707 break;
708 case s_D2:
709 pointMagStart = HepGeom::Point3D<double>(m_Config.pointD2Start[(beam-1)*3],m_Config.pointD2Start[(beam-1)*3+1],m_Config.pointD2Start[(beam-1)*3+2]);
710 pointMagEnd = HepGeom::Point3D<double>(m_Config.pointD2End[(beam-1)*3],m_Config.pointD2End[(beam-1)*3+1],m_Config.pointD2End[(beam-1)*3+2]);
711 rotZ = m_Config.fD2RotZ[beam-1];
712 break;
713 }
714}
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:63
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
Point
ChargedTracksWeightFilter::Spline::Point Point
Definition ChargedTracksWeightFilter.cxx:11
length
double length(const pvec &v)
Definition FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
PathResolverFindDataFile
std::string PathResolverFindDataFile(const std::string &logical_file_name)
Definition PathResolver.cxx:320
angle
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
Definition TRTDetectorFactory_Full.cxx:71
y
#define y
xyz
#define xyz
x
#define x
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15
MagField::ForwardRegionFieldSvc::loadDataFileNLines
std::vector< std::vector< std::string > > loadDataFileNLines(const char *fileName, int cols, int N, std::string &header)
Definition ForwardRegionFieldSvc.cxx:423
MagField::ForwardRegionFieldSvc::m_magIndMQXB
double m_magIndMQXB[s_rowsMQXB][s_colsMQXB][2]
Definition ForwardRegionFieldSvc.h:118
MagField::ForwardRegionFieldSvc::getMagNumFromName
int getMagNumFromName(const std::string &name) const
Definition ForwardRegionFieldSvc.cxx:580
MagField::ForwardRegionFieldSvc::m_magnets
std::vector< std::vector< std::string > > m_magnets
Definition ForwardRegionFieldSvc.h:75
MagField::ForwardRegionFieldSvc::getField
virtual void getField(const double *xyz, double *bxyz, double *deriv=nullptr) const override final
IMagFieldSvc interface methods.
Definition ForwardRegionFieldSvc.cxx:123
MagField::ForwardRegionFieldSvc::getMag
double getMag(int magnet, int beam) const
Definition ForwardRegionFieldSvc.cxx:398
MagField::ForwardRegionFieldSvc::m_magIndMQXA
double m_magIndMQXA[s_rowsMQXA][s_colsMQXA][2]
Definition ForwardRegionFieldSvc.h:117
MagField::ForwardRegionFieldSvc::ForwardRegionFieldSvc
ForwardRegionFieldSvc(const std::string &name, ISvcLocator *svc)
Constructor with parameters.
Definition ForwardRegionFieldSvc.cxx:37
MagField::ForwardRegionFieldSvc::writeOutTwiss
void writeOutTwiss(const std::vector< std::vector< std::string > > &loadedTwissFile, int beam, const std::string &header)
Definition ForwardRegionFieldSvc.cxx:519
MagField::ForwardRegionFieldSvc::handle
void handle(const Incident &runIncident) override final
IIncidentListener interface methods.
Definition ForwardRegionFieldSvc.cxx:108
MagField::ForwardRegionFieldSvc::m_properties
ToolHandle< IForwardRegionProperties > m_properties
Definition ForwardRegionFieldSvc.h:132
MagField::ForwardRegionFieldSvc::loadDataFile
std::vector< std::vector< std::string > > loadDataFile(const char *fileName, int cols)
Definition ForwardRegionFieldSvc.cxx:507
MagField::ForwardRegionFieldSvc::FieldValue
G4ThreeVector FieldValue(G4ThreeVector Point) const
Definition ForwardRegionFieldSvc.cxx:140
MagField::ForwardRegionFieldSvc::InitMagDataFromTwiss
void InitMagDataFromTwiss(const std::vector< std::vector< std::string > > &loadedTwissFile, int beam, double momentum)
Definition ForwardRegionFieldSvc.cxx:250
MagField::ForwardRegionFieldSvc::kLToB
double kLToB(double kL, double lenght, double momentum)
Definition ForwardRegionFieldSvc.cxx:514
MagField::ForwardRegionFieldSvc::initialize
StatusCode initialize() override final
AthService interface methods.
Definition ForwardRegionFieldSvc.cxx:95
MagField::ForwardRegionFieldSvc::getMagnetTransformParams
void getMagnetTransformParams(int beam, int magnet, G4ThreeVector Point, HepGeom::Point3D< double > &pointMagStart, HepGeom::Point3D< double > &pointMagEnd, double &rotZ) const
Definition ForwardRegionFieldSvc.cxx:617
MagField::ForwardRegionFieldSvc::getMagInd
G4ThreeVector getMagInd(G4ThreeVector Point, int q, int beam) const
Non-inherited public methods FIXME - add new interface?
Definition ForwardRegionFieldSvc.cxx:324
MagField::ForwardRegionFieldSvc::getKick
G4ThreeVector getKick(int beam) const
Definition ForwardRegionFieldSvc.cxx:404
MagField::ForwardRegionFieldSvc::getFieldZR
virtual void getFieldZR(const double *xyz, double *bxyz, double *deriv=nullptr) const override final
get B field value on the z-r plane at given position
Definition ForwardRegionFieldSvc.cxx:135
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition getMessageSvc.cxx:20
file
TFile * file
Definition tile_monitor.h:29
Generated on for ATLAS Offline Software by doxygen 1.14.0