BmnRoot
Loading...
Searching...
No Matches
BmnSP41FieldMapCreator.cxx
Go to the documentation of this file.
1#include "BmnSP41FieldMapCreator.h"
2#include "BmnNewFieldMap.h"
3
4#include "TFitResult.h"
5#include "TCanvas.h"
6#include "TStyle.h"
7#include "TLorentzVector.h"
8
9using namespace std;
10
11fitParMonitor::fitParMonitor(TString out)
12: C(nullptr),
13 M(nullptr),
14 S(nullptr),
15 f(nullptr)
16{
17 const Int_t nParts = 2; // Left part, right part of extrapolation to be done ...
18 const Int_t nComp = 3; // Left part, right part of extrapolation to be done ...
19
20 C = new TH1F**[nParts];
21 M = new TH1F**[nParts];
22 S = new TH1F**[nParts];
23
24 f = new TFile(Form("fitParMonitor_%s.root", out.Data()), "recreate");
25
26 for (Int_t iPart = 0; iPart < nParts; iPart++) {
27 C[iPart] = new TH1F*[nComp];
28 M[iPart] = new TH1F*[nComp];
29 S[iPart] = new TH1F*[nComp];
30 for (Int_t iComp = 0; iComp < nComp; iComp++) {
31 C[iPart][iComp] = new TH1F(Form("Constant_%d_%d", iPart, iComp), "Constant", 100, 0., 0.);
32 M[iPart][iComp] = new TH1F(Form("Mean_%d_%d", iPart, iComp), "Mean", 100, 0., 0.);
33 S[iPart][iComp] = new TH1F(Form("Sigma_%d_%d", iPart, iComp), "Sigma", 100, 0., 0.);
34 }
35 }
36}
37
38fitParMonitor::~fitParMonitor()
39{
40 for (Int_t iPart = 0; iPart < 2; iPart++) {
41 for (Int_t iComp = 0; iComp < 3; iComp++) {
42 C[iPart][iComp]->Write();
43 M[iPart][iComp]->Write();
44 S[iPart][iComp]->Write();
45 }
46 }
47
48 if (f)
49 delete f;
50}
51
52BmnSP41FieldMapCreator::~BmnSP41FieldMapCreator() {
53 for (Int_t iPart = 0; iPart < 2; iPart++) {
54 for (Int_t iComponent = 0; iComponent < 3; iComponent++) {
55 delete constant[iPart][iComponent];
56 delete mean[iPart][iComponent];
57 delete sigma[iPart][iComponent];
58 }
59 delete [] constant[iPart];
60 delete [] mean[iPart];
61 delete [] sigma[iPart];
62 }
63 delete [] constant;
64 delete [] mean;
65 delete [] sigma;
66
67
68 if (fMap)
69 delete fMap;
70
71 if (fConditions)
72 delete fConditions;
73
74 if (fMonitor)
75 delete fMonitor;
76}
77
78BmnSP41FieldMapCreator::BmnSP41FieldMapCreator(TString in, TString out)
79: fDebug(kFALSE),
80 nameOutput(out),
81 dimGrid(200),
82 dimNhist(6),
83 factor(10.)
84{
85 // Read initial field map
86 fMap = new BmnNewFieldMap(in);
87 fMap->SetScale(1.);
88 fMap->Init();
89
90 // Extrapolation conditions ...
91 fConditions = new extrapolationConditions();
92
93 // Fit par monitor ...
94 fMonitor = new fitParMonitor(out);
95
96 WinXmin = fMap->GetXmin();
97 WinXmax = fMap->GetXmax();
98
99 WinYmin = fMap->GetYmin();
100 WinYmax = fMap->GetYmax();
101
102 WinZmin = fMap->GetZmin();
103 WinZmax = fMap->GetZmax();
104
105 winXstep = (WinXmax - WinXmin) / (dimNhist - 1);
106 winYstep = (WinYmax - WinYmin) / (dimNhist - 1);
107 winZstep = (WinZmax - WinZmin) / (dimNhist - 1);
108
109 const Int_t nComponents = 3;
110 const Int_t nParts = 2;
111
112 constant = new TGraph2D**[nParts];
113 mean = new TGraph2D**[nParts];
114 sigma = new TGraph2D**[nParts];
115 for (Int_t iPart = 0; iPart < nParts; iPart++) {
116 constant[iPart] = new TGraph2D*[nComponents];
117 mean[iPart] = new TGraph2D*[nComponents];
118 sigma[iPart] = new TGraph2D*[nComponents];
119 for (Int_t iComponent = 0; iComponent < nComponents; iComponent++) {
120 constant[iPart][iComponent] = new TGraph2D();
121 mean[iPart][iComponent] = new TGraph2D();
122 sigma[iPart][iComponent] = new TGraph2D();
123 }
124 }
125}
126
127void BmnSP41FieldMapCreator::Extrapolation() {
128 TString axis = fConditions->direction;
129
130 // Setting predefined tail ranges (left and right) ...
131 // They are considered as subjects to change if necessary ...
132 if (axis.Contains("X")) {
133 fConditions->tailRangeLeft = make_pair(-100., -80.);
134 fConditions->tailRangeRight = make_pair(80., 100.);
135 } else if (axis.Contains("Z")) {
136 fConditions->tailRangeLeft = make_pair(-28., 0.);
137 fConditions->tailRangeRight = make_pair(340., 360.);
138 }
139
140 const Int_t nParts = 2;
141 const Int_t nComponents = 3;
142
143 if (axis.Contains("X")) {
144 for (Int_t iPart = 0; iPart < nParts; iPart++)
145 for (Int_t iComp = 0; iComp < nComponents; iComp++)
146 DoExtrapolationAlongX(iPart, iComp);
147 } else if (axis.Contains("Z")) {
148 for (Int_t iPart = 0; iPart < nParts; iPart++)
149 for (Int_t iComp = 0; iComp < nComponents; iComp++)
150 DoExtrapolationAlongZ(iPart, iComp);
151 } else {
152 cout << "Not supported direction!!! Exiting ..." << endl;
153 throw;
154 }
155
156 CreateExtrapolatedMap();
157}
158
159void BmnSP41FieldMapCreator::DoExtrapolationAlongX(Int_t iPart, Int_t iComp) {
160 TGraph*** bField = new TGraph**[dimGrid];
161 TF1*** fit_func = new TF1**[dimNhist];
162 TString comp = (iComp == 0) ? "x" : (iComp == 1) ? "y" : (iComp == 2) ? "z" : "";
163
164 Double_t Ymax = fMap->GetYmax();
165 Double_t Ymin = fMap->GetYmin();
166
167 Double_t Zmax = fMap->GetZmax();
168 Double_t Zmin = fMap->GetZmin();
169
170 // Define number of slices on Y(Z)-dimension of the grid used
171 Int_t nSlicesY = dimGrid - 1;
172 Int_t nSlicesZ = dimGrid - 1;
173
174 // Define step along X(Y, Z)-axises of the used grid
175 Double_t dz = (Zmax - Zmin) / nSlicesZ;
176 Double_t dy = (Ymax - Ymin) / nSlicesY;
177
178 Int_t counter = 0;
179
180 for (Int_t i = 0; i < dimGrid; i++) {
181 if (fDebug)
182 if (i % 10 == 0)
183 cout << "Processing node " << i << endl;
184 bField[i] = new TGraph*[dimGrid];
185 Double_t ySlice = Ymin + i * dy;
186
187 for (Int_t j = 0; j < dimGrid; j++) {
188 Float_t zSlice = Zmin + j * dz;
189 bField[i][j] = new TGraph();
190
191 for (Int_t k = 0; k < fMap->GetNx(); k++) {
192 Double_t x = fMap->GetXmin() + k * fMap->GetXstep();
193 bField[i][j]->SetPoint(k, x, (comp == "x") ? fMap->GetBx(x, ySlice, zSlice) :
194 (comp == "y") ? fMap->GetBy(x, ySlice, zSlice) :
195 (comp == "z") ? fMap->GetBz(x, ySlice, zSlice) : 0.0);
196 }
197
198 vector <Double_t> fit_par;
199
200 TF1* f = nullptr;
201 if (iPart == 0) {
202 f = new TF1("f", "gaus", fConditions->tailRangeLeft.first, fConditions->tailRangeLeft.second);
203 f->SetParameter(2, 50.);
204 fit_par = bField[i][j]->Fit(f, "SQ", "", f->GetXmin(), f->GetXmax())->Parameters();
205 } else {
206 f = new TF1("f", "gaus", fConditions->tailRangeRight.first, fConditions->tailRangeRight.second);
207 f->SetParameter(2, 100.);
208 fit_par = bField[i][j]->Fit(f, "SQ", "", f->GetXmin(), f->GetXmax())->Parameters();
209 }
210
211 fMonitor->C[iPart][iComp]->Fill(fit_par[0]);
212 fMonitor->M[iPart][iComp]->Fill(fit_par[1]);
213 fMonitor->S[iPart][iComp]->Fill(fit_par[2]);
214
215 constant[iPart][iComp]->SetPoint(counter, ySlice, zSlice, fit_par[0]);
216 mean[iPart][iComp]->SetPoint(counter, ySlice, zSlice, fit_par[1]);
217 sigma[iPart][iComp]->SetPoint(counter, ySlice, zSlice, fit_par[2]);
218
219 if (f)
220 delete f;
221
222 counter++;
223 }
224 }
225
226 for (Int_t i = 0; i < dimNhist; i++) {
227 fit_func[i] = new TF1*[dimNhist];
228 Double_t y = WinYmin + i * winYstep;
229
230 for (Int_t j = 0; j < dimNhist; j++) {
231 Double_t z = WinZmin + j * winZstep;
232 fit_func[i][j] = new TF1("func", "[0]*exp(-0.5*((x-[1])/[2])**2)",
233 (iPart == 1) ? fConditions->tailRangeRight.second : fConditions->min, (iPart == 1) ? fConditions->max : fConditions->tailRangeLeft.first);
234 fit_func[i][j]->SetParameter(0, constant[iPart][iComp]->Interpolate(y, z));
235 fit_func[i][j]->SetParameter(1, mean[iPart][iComp]->Interpolate(y, z));
236 fit_func[i][j]->SetParameter(2, sigma[iPart][iComp]->Interpolate(y, z));
237 }
238 }
239
240 DrawProfileExtrapX(dimNhist, bField, fit_func, "B" + comp + ".png", comp, Ymin, Zmin, dy, dz);
241
242 for (Int_t i = 0; i < dimNhist; i++) {
243 for (Int_t j = 0; j < dimNhist; j++)
244 delete fit_func[i][j];
245 delete [] fit_func[i];
246 }
247 delete [] fit_func;
248
249 for (Int_t i = 0; i < dimGrid; i++) {
250 for (Int_t j = 0; j < dimGrid; j++)
251 delete bField[i][j];
252 delete [] bField[i];
253 }
254 delete [] bField;
255}
256
257void BmnSP41FieldMapCreator::DoExtrapolationAlongZ(Int_t iPart, Int_t iComp) {
258 TGraph*** bField = new TGraph**[dimGrid];
259 TF1*** fit_func = new TF1**[dimNhist];
260 TString comp = (iComp == 0) ? "x" : (iComp == 1) ? "y" : (iComp == 2) ? "z" : "";
261
262 Double_t Xmax = fMap->GetXmax();
263 Double_t Xmin = fMap->GetXmin();
264
265 Double_t Ymax = fMap->GetYmax();
266 Double_t Ymin = fMap->GetYmin();
267
268 // Define number of slices on X(Y)-dimension of the grid used
269 Int_t nSlicesX = dimGrid - 1;
270 Int_t nSlicesY = dimGrid - 1;
271
272 // Define step along X(Y, Z)-axises of the used grid
273 Double_t dx = (Xmax - Xmin) / nSlicesX;
274 Double_t dy = (Ymax - Ymin) / nSlicesY;
275
276 Int_t counter = 0;
277
278 for (Int_t i = 0; i < dimGrid; i++) {
279 if (fDebug)
280 if (i % 10 == 0)
281 cout << "Processing node " << i << endl;
282 bField[i] = new TGraph*[dimGrid];
283 Double_t xSlice = Xmin + i * dx;
284
285 for (Int_t j = 0; j < dimGrid; j++) {
286 Float_t ySlice = Ymin + j * dy;
287 bField[i][j] = new TGraph();
288
289 for (Int_t k = 0; k < fMap->GetNz(); k++) {
290 Double_t z = fMap->GetZmin() + k * fMap->GetZstep();
291 bField[i][j]->SetPoint(k, z, (comp == "x") ? fMap->GetBx(xSlice, ySlice, z) :
292 (comp == "y") ? fMap->GetBy(xSlice, ySlice, z) :
293 (comp == "z") ? fMap->GetBz(xSlice, ySlice, z) : 0.0);
294 }
295
296 vector <Double_t> fit_par;
297
298 TF1* f = nullptr;
299 if (iPart == 0) {
300 f = new TF1("f", "gaus", fConditions->tailRangeLeft.first, fConditions->tailRangeLeft.second);
301 f->SetParameter(2, 100.);
302 fit_par = bField[i][j]->Fit(f, "SQ", "G", f->GetXmin(), f->GetXmax())->Parameters();
303 } else {
304 f = new TF1("f", "gaus", fConditions->tailRangeRight.first, fConditions->tailRangeRight.second);
305 f->SetParameter(2, 120.);
306 fit_par = bField[i][j]->Fit(f, "SQ", "G", f->GetXmin(), f->GetXmax())->Parameters();
307 }
308
309 fMonitor->C[iPart][iComp]->Fill(fit_par[0]);
310 fMonitor->M[iPart][iComp]->Fill(fit_par[1]);
311 fMonitor->S[iPart][iComp]->Fill(fit_par[2]);
312
313 constant[iPart][iComp]->SetPoint(counter, xSlice, ySlice, fit_par[0]);
314 mean[iPart][iComp]->SetPoint(counter, xSlice, ySlice, fit_par[1]);
315 sigma[iPart][iComp]->SetPoint(counter, xSlice, ySlice, fit_par[2]);
316
317 counter++;
318
319 if (f)
320 delete f;
321 }
322 }
323
324 for (Int_t i = 0; i < dimNhist; i++) {
325 fit_func[i] = new TF1*[dimNhist];
326 Double_t x = WinXmin + i * winXstep;
327
328 for (Int_t j = 0; j < dimNhist; j++) {
329 Double_t y = WinYmin + j * winYstep;
330 fit_func[i][j] = new TF1("func", "[0]*exp(-0.5*((x-[1])/[2])**2)", (iPart == 1) ? fConditions->tailRangeRight.second : fConditions->min,
331 (iPart == 1) ? fConditions->max : fConditions->tailRangeLeft.first);
332 fit_func[i][j]->SetParameter(0, constant[iPart][iComp]->Interpolate(x, y));
333 fit_func[i][j]->SetParameter(1, mean[iPart][iComp]->Interpolate(x, y));
334 fit_func[i][j]->SetParameter(2, sigma[iPart][iComp]->Interpolate(x, y));
335 }
336 }
337
338 DrawProfileExtrapZ(dimNhist, bField, fit_func, "B" + comp + ".pdf", comp, Xmin, Ymin, dx, dy);
339
340 for (Int_t i = 0; i < dimNhist; i++) {
341 for (Int_t j = 0; j < dimNhist; j++)
342 delete fit_func[i][j];
343 delete [] fit_func[i];
344 }
345 delete [] fit_func;
346
347 for (Int_t i = 0; i < dimGrid; i++) {
348 for (Int_t j = 0; j < dimGrid; j++)
349 delete bField[i][j];
350 delete [] bField[i];
351 }
352 delete [] bField;
353}
354
355void BmnSP41FieldMapCreator::DrawProfileExtrapZ(Int_t dim, TGraph*** gr, TF1*** f, TString name, TString comp, Double_t Xmin, Double_t Ymin, Double_t dx, Double_t dy) {
356 TCanvas* c = new TCanvas("1", "1", 1200, 800);
357 c->Divide(dim, dim);
358 Int_t counter = 1;
359 for (Int_t i = 0; i < dim; i++) {
360 Double_t X = WinXmin + i * winXstep;
361 Double_t sliceX = (X - Xmin) / dx;
362 Int_t shiftI = Int_t(sliceX);
363 for (Int_t j = 0; j < dim; j++) {
364 Double_t Y = WinYmin + j * winYstep;
365 Double_t sliceY = (Y - Ymin) / dy;
366 Int_t shiftJ = Int_t(sliceY);
367 c->cd(counter);
368 gr[shiftI][shiftJ]->Draw();
369 gr[shiftI][shiftJ]->GetXaxis()->SetLimits(fConditions->min, fConditions->max);
370 gr[shiftI][shiftJ]->SetTitle("B" + comp + TString(Form(", x = %G [cm]", X)) + TString(Form(", y = %G [cm]", Y)));
371 gr[shiftI][shiftJ]->GetXaxis()->SetTitleSize(0.1);
372 gr[shiftI][shiftJ]->GetXaxis()->SetLabelSize(0.08);
373 gr[shiftI][shiftJ]->GetYaxis()->SetTitleSize(0.1);
374 gr[shiftI][shiftJ]->GetYaxis()->SetLabelSize(0.1);
375 gr[shiftI][shiftJ]->GetXaxis()->SetTitle("Z, cm");
376 gr[shiftI][shiftJ]->GetYaxis()->SetTitle("B" + comp + ", kG");
377 counter++;
378 f[i][j]->Draw("same");
379 }
380 }
381 c->SaveAs(name);
382 delete c;
383}
384
385void BmnSP41FieldMapCreator::DrawProfileExtrapX(Int_t dim, TGraph*** gr, TF1*** f, TString name, TString comp, Double_t Ymin, Double_t Zmin, Double_t dy, Double_t dz) {
386 gStyle->SetTitleFontSize(0.1);
387 TCanvas* c = new TCanvas("1", "1", 4000, 4000);
388 c->Divide(dim, dim);
389 Int_t counter = 1;
390 for (Int_t i = 0; i < dim; i++) {
391 Double_t Y = WinYmin + i * winYstep;
392 Double_t sliceY = (Y - Ymin) / dy;
393 Int_t shiftI = Int_t(sliceY);
394 for (Int_t j = 0; j < dim; j++) {
395 Double_t Z = WinZmin + j * winZstep;
396 Double_t sliceZ = (Z - Zmin) / dz;
397 Int_t shiftJ = Int_t(sliceZ);
398 c->cd(counter);
399 gr[shiftI][shiftJ]->Draw();
400 gr[shiftI][shiftJ]->GetXaxis()->SetLimits(fConditions->min, fConditions->max);
401 gr[shiftI][shiftJ]->SetLineWidth(2);
402 gr[shiftI][shiftJ]->GetXaxis()->SetNdivisions(5);
403 gr[shiftI][shiftJ]->SetTitle("B" + comp + TString(Form(", y = %G [cm]", Y)) + TString(Form(", z = %G [cm]", Z)));
404 gr[shiftI][shiftJ]->GetXaxis()->SetTitleSize(0.1);
405 gr[shiftI][shiftJ]->GetXaxis()->SetLabelSize(0.08);
406 gr[shiftI][shiftJ]->GetYaxis()->SetTitleSize(0.1);
407 gr[shiftI][shiftJ]->GetYaxis()->SetLabelSize(0.1);
408 gr[shiftI][shiftJ]->GetXaxis()->SetTitle("X, cm");
409 gr[shiftI][shiftJ]->GetYaxis()->SetTitle("B" + comp + ", kG");
410 counter++;
411 f[i][j]->Draw("same");
412 f[i][j]->SetLineWidth(2);
413 }
414 }
415 c->SaveAs(name);
416 delete c;
417}
418
419void BmnSP41FieldMapCreator::CreateExtrapolatedMap() {
420
421 Double_t Xmin = fMap->GetXmin();
422 Double_t Ymin = fMap->GetYmin();
423 Double_t Zmin = fMap->GetZmin();
424
425 Double_t Xmax = fMap->GetXmax();
426 Double_t Ymax = fMap->GetYmax();
427 Double_t Zmax = fMap->GetZmax();
428
429 Double_t dx = fMap->GetXstep();
430 Double_t dy = fMap->GetYstep();
431 Double_t dz = fMap->GetZstep();
432
433 Int_t Nx = fMap->GetNx();
434 Int_t Ny = fMap->GetNy();
435 Int_t Nz = fMap->GetNz();
436 // Int_t Nz = (Zmax - Zmin) / dz + 1;
437
438 // Print(fMap);
439
440 Double_t x = -1., y = -1., z = -1.;
441
442 ofstream mapFile(nameOutput);
443
444 mapFile.precision(6);
445 mapFile << showpoint;
446 mapFile << "nosym" << endl;
447
448 TString axis = fConditions->direction;
449
450 // Filling parameters of the extrapolated map ...
451 if (axis.Contains("X")) {
452 Xmin = fConditions->min;
453 Xmax = fConditions->max;
454 Nx = Int_t((Xmax - Xmin) / 2.) + 1;
455 mapFile << Xmin << " " << Xmax << " " << Nx << endl;
456 } else
457 mapFile << Xmin << " " << Xmax << " " << Nx << endl;
458
459 mapFile << Ymin << " " << Ymax << " " << Ny << endl;
460
461 if (axis.Contains("Z")) {
462 Zmin = fConditions->min;
463 Zmax = fConditions->max;
464 Nz = Int_t((Zmax - Zmin) / 2.) + 1;
465 mapFile << Zmin << " " << Zmax << " " << Nz << endl;
466 } else
467 mapFile << Zmin << " " << Zmax << " " << Nz << endl;
468
469 for (Int_t ix = 0; ix < Nx; ix++) {
470 x = Xmin + dx * ix;
471 for (Int_t iy = 0; iy < Ny; iy++) {
472 y = Ymin + dy * iy;
473 if (fDebug)
474 cout << x << " " << y << endl;
475 for (Int_t iz = 0; iz < Nz; iz++) {
476 //Int_t ind = ix * Ny * Nz + iy * Nz + iz;
477 z = Zmin + dz * iz;
478
479 Double_t bx = -1., by = -1., bz = -1.;
480
481 Bool_t isInMap = kTRUE;
482
483 if (axis.Contains("X")) {
484 if (x < fMap->GetXmin()) {
485 isInMap = kFALSE;
486 bx = FieldExtrapolate(x, y, z, constant[0][0], mean[0][0], sigma[0][0]);
487 by = FieldExtrapolate(x, y, z, constant[0][1], mean[0][1], sigma[0][1]);
488 bz = FieldExtrapolate(x, y, z, constant[0][2], mean[0][2], sigma[0][2]);
489 } else if (x > fMap->GetXmax() - 1.) {
490 isInMap = kFALSE;
491 bx = FieldExtrapolate(x, y, z, constant[1][0], mean[1][0], sigma[1][0]);
492 by = FieldExtrapolate(x, y, z, constant[1][1], mean[1][1], sigma[1][1]);
493 bz = FieldExtrapolate(x, y, z, constant[1][2], mean[1][2], sigma[1][2]);
494 }
495 } else if (axis.Contains("Z")) {
496 if (z < fMap->GetZmin()) {
497 isInMap = kFALSE;
498 bx = FieldExtrapolate(x, y, z, constant[0][0], mean[0][0], sigma[0][0]);
499 by = FieldExtrapolate(x, y, z, constant[0][1], mean[0][1], sigma[0][1]);
500 bz = FieldExtrapolate(x, y, z, constant[0][2], mean[0][2], sigma[0][2]);
501 } else if (z > fMap->GetZmax() - 1.) {
502 isInMap = kFALSE;
503 bx = FieldExtrapolate(x, y, z, constant[1][0], mean[1][0], sigma[1][0]);
504 by = FieldExtrapolate(x, y, z, constant[1][1], mean[1][1], sigma[1][1]);
505 bz = FieldExtrapolate(x, y, z, constant[1][2], mean[1][2], sigma[1][2]);
506 }
507 }
508
509 if (isInMap) {
510 bx = fMap->GetBx(x, y, z);
511 by = fMap->GetBy(x, y, z);
512 bz = fMap->GetBz(x, y, z);
513 }
514
515 // Checking obtained extrapolated values ...
516 mapFile << checkedFieldValue(bx, by, bz).at(0) / factor << " " << checkedFieldValue(bx, by, bz).at(1) / factor << " " <<
517 checkedFieldValue(bx, by, bz).at(2) / factor << endl;
518 } // z-Loop
519 } // y-Loop
520 } // x-Loop
521}
522
523void BmnSP41FieldMapCreator::SmoothMap(TClonesArray* sValues) {
524
525 vector <TLorentzVector> vSmoo;
526
527 for (Int_t iSmooth = 0; iSmooth < sValues->GetEntriesFast(); iSmooth++) {
528 smoothedValues* sValue = (smoothedValues*) sValues->UncheckedAt(iSmooth);
529
530 TLorentzVector tmp(sValue->X, sValue->Y, sValue->Z, sValue->fieldValue);
531 vSmoo.push_back(tmp);
532 }
533
534 Double_t Xmin = fMap->GetXmin();
535 Double_t Ymin = fMap->GetYmin();
536 Double_t Zmin = fMap->GetZmin();
537 Double_t Zmax = fMap->GetZmax();
538
539 Double_t dx = fMap->GetXstep();
540 Double_t dy = fMap->GetYstep();
541 Double_t dz = fMap->GetZstep();
542
543 Int_t Nx = fMap->GetNx();
544 Int_t Ny = fMap->GetNy();
545 Int_t Nz = fMap->GetNz();
546
547 Print(fMap);
548
549 Double_t x, y, z, bx, by, bz;
550
551 UInt_t index = 0;
552 ofstream mapFile(nameOutput);
553
554 mapFile.precision(6);
555 mapFile << showpoint;
556 mapFile << "nosym" << endl;
557 mapFile << Xmin << " " << -Xmin << " " << Nx << endl;
558 mapFile << Ymin << " " << -Ymin << " " << Ny << endl;
559 mapFile << Zmin << " " << Zmax << " " << Nz << endl;
560
561 for (Int_t ix = 0; ix < Nx; ix++) {
562 x = Xmin + dx * ix;
563
564 for (Int_t iy = 0; iy < Ny; iy++) {
565 y = Ymin + dy * iy;
566
567 for (Int_t iz = 0; iz < Nz; iz++) {
568 index = ix * Ny * Nz + iy * Nz + iz;
569 z = Zmin + dz * iz;
570
571 by = (fMap->GetBy(x, y, z) / factor);
572 for (auto it = vSmoo.begin(); it != vSmoo.end(); it++) {
573
574 TLorentzVector v = *it;
575
576 if (v.X() != x || v.Y() != y || v.Z() != z)
577 continue;
578
579 by = v.T() / factor;
580
581 vSmoo.erase(it--);
582 break;
583 }
584
585 bx = (Abs(fMap->GetBx(x, y, z)) / factor > 10e-8) ? (fMap->GetBx(x, y, z) / factor) : 0.0;
586 bz = (Abs(fMap->GetBz(x, y, z)) / factor > 10e-8) ? (fMap->GetBz(x, y, z) / factor) : 0.0;
587 // mapFile << x << " " << y << " " << z << " " << bx << " " << by << " " << bz << endl;
588 mapFile << bx << " " << by << " " << bz << endl;
589
590 if (index % 1000 == 0)
591 cout << "nNodes processed# " << index << endl;
592
593 } // z-Loop
594 } // y-Loop
595 } // x-Loop
596}
597
598Double_t BmnSP41FieldMapCreator::FieldExtrapolate(Double_t x, Double_t y, Double_t z, TGraph2D* c, TGraph2D* m, TGraph2D* s) {
599 vector <Double_t> varToBePassed;
600
601 if (fConditions->direction.Contains("X")) {
602 varToBePassed.push_back(y);
603 varToBePassed.push_back(z);
604 varToBePassed.push_back(x);
605 } else if (fConditions->direction.Contains("Z")) {
606 varToBePassed.push_back(x);
607 varToBePassed.push_back(y);
608 varToBePassed.push_back(z);
609 }
610
611 Double_t Mean = m->Interpolate(varToBePassed[0], varToBePassed[1]);
612 Double_t Sigma = s->Interpolate(varToBePassed[0], varToBePassed[1]);
613
614 return c->Interpolate(varToBePassed[0], varToBePassed[1]) * TMath::Exp(-0.5 * (varToBePassed[2] - Mean) * (varToBePassed[2] - Mean) / Sigma / Sigma);
615}
616
617void BmnSP41FieldMapCreator::Print(BmnFieldMap* magField) {
618 cout << "Nx = " << magField->GetNx();
619 cout << "; Ny = " << magField->GetNy();
620 cout << "; Nz = " << magField->GetNz();
621
622 Double_t shift = 1e-6;
623
624 cout << "; Xmin = " << magField->GetXmin() - shift;
625 cout << "; dx = " << magField->GetXstep();
626 cout << "; Ymin = " << magField->GetYmin() - shift;
627 cout << "; dy = " << magField->GetYstep();
628 cout << "; Zmin = " << magField->GetZmin() - shift;
629 cout << "; dz = " << magField->GetZstep() << endl;
630}
Sigma
Double_t Sigma(vector< Double_t > dist, vector< Double_t > w)
Definition BmnMath.cxx:890
v
__m128 v
Definition P4_F32vec4.h:1
i
int i
Definition P4_F32vec4.h:22
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:30
m
__m128 m
Definition P4_F32vec4.h:27
f
float f
Definition P4_F32vec4.h:21
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:31
BmnFieldMap::GetNz
Int_t GetNz() const
Definition BmnFieldMap.h:82
BmnFieldMap::GetBz
virtual Double_t GetBz(Double_t x, Double_t y, Double_t z)=0
BmnFieldMap::GetYstep
Double_t GetYstep() const
Definition BmnFieldMap.h:77
BmnFieldMap::GetZstep
Double_t GetZstep() const
Definition BmnFieldMap.h:78
BmnFieldMap::GetXstep
Double_t GetXstep() const
Definition BmnFieldMap.h:76
BmnFieldMap::GetBy
virtual Double_t GetBy(Double_t x, Double_t y, Double_t z)=0
BmnFieldMap::GetZmin
Double_t GetZmin() const
Definition BmnFieldMap.h:70
BmnFieldMap::GetYmin
Double_t GetYmin() const
Definition BmnFieldMap.h:69
BmnFieldMap::SetScale
void SetScale(Double_t factor)
Definition BmnFieldMap.h:55
BmnFieldMap::GetBx
virtual Double_t GetBx(Double_t x, Double_t y, Double_t z)=0
BmnFieldMap::GetXmax
Double_t GetXmax() const
Definition BmnFieldMap.h:72
BmnFieldMap::GetNy
Int_t GetNy() const
Definition BmnFieldMap.h:81
BmnFieldMap::GetNx
Int_t GetNx() const
Definition BmnFieldMap.h:80
BmnFieldMap::Init
virtual void Init()
Definition BmnFieldMap.cxx:93
BmnFieldMap::GetXmin
Double_t GetXmin() const
Definition BmnFieldMap.h:68
BmnFieldMap::GetZmax
Double_t GetZmax() const
Definition BmnFieldMap.h:74
BmnFieldMap::GetYmax
Double_t GetYmax() const
Definition BmnFieldMap.h:73
extrapolationConditions::tailRangeRight
std::pair< Double_t, Double_t > tailRangeRight
Definition BmnSP41FieldMapCreator.h:54
extrapolationConditions::tailRangeLeft
std::pair< Double_t, Double_t > tailRangeLeft
Definition BmnSP41FieldMapCreator.h:53
dim
const UInt_t dim
Definition MpdVHLLEGenerator.h:21
check_success.s
s
Definition check_success.py:47
std
STL namespace.