BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnSP41FieldMapCreator.cxx
Go to the documentation of this file.
1#include "BmnSP41FieldMapCreator.h"
2
3#include "BmnNewFieldMap.h"
4#include "TCanvas.h"
5#include "TFitResult.h"
6#include "TLorentzVector.h"
7#include "TStyle.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{
54 for (Int_t iPart = 0; iPart < 2; iPart++) {
55 for (Int_t iComponent = 0; iComponent < 3; iComponent++) {
56 delete constant[iPart][iComponent];
57 delete mean[iPart][iComponent];
58 delete sigma[iPart][iComponent];
59 }
60 delete[] constant[iPart];
61 delete[] mean[iPart];
62 delete[] sigma[iPart];
63 }
64 delete[] constant;
65 delete[] mean;
66 delete[] sigma;
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{
129 TString axis = fConditions->direction;
130
131 // Setting predefined tail ranges (left and right) ...
132 // They are considered as subjects to change if necessary ...
133 if (axis.Contains("X")) {
134 fConditions->tailRangeLeft = make_pair(-100., -80.);
135 fConditions->tailRangeRight = make_pair(80., 100.);
136 } else if (axis.Contains("Z")) {
137 fConditions->tailRangeLeft = make_pair(-28., 0.);
138 fConditions->tailRangeRight = make_pair(340., 360.);
139 }
140
141 const Int_t nParts = 2;
142 const Int_t nComponents = 3;
143
144 if (axis.Contains("X")) {
145 for (Int_t iPart = 0; iPart < nParts; iPart++)
146 for (Int_t iComp = 0; iComp < nComponents; iComp++)
147 DoExtrapolationAlongX(iPart, iComp);
148 } else if (axis.Contains("Z")) {
149 for (Int_t iPart = 0; iPart < nParts; iPart++)
150 for (Int_t iComp = 0; iComp < nComponents; iComp++)
151 DoExtrapolationAlongZ(iPart, iComp);
152 } else {
153 cout << "Not supported direction!!! Exiting ..." << endl;
154 throw;
155 }
156
157 CreateExtrapolatedMap();
158}
159
160void BmnSP41FieldMapCreator::DoExtrapolationAlongX(Int_t iPart, Int_t iComp)
161{
162 TGraph*** bField = new TGraph**[dimGrid];
163 TF1*** fit_func = new TF1**[dimNhist];
164 TString comp = (iComp == 0) ? "x" : (iComp == 1) ? "y" : (iComp == 2) ? "z" : "";
165
166 Double_t Ymax = fMap->GetYmax();
167 Double_t Ymin = fMap->GetYmin();
168
169 Double_t Zmax = fMap->GetZmax();
170 Double_t Zmin = fMap->GetZmin();
171
172 // Define number of slices on Y(Z)-dimension of the grid used
173 Int_t nSlicesY = dimGrid - 1;
174 Int_t nSlicesZ = dimGrid - 1;
175
176 // Define step along X(Y, Z)-axises of the used grid
177 Double_t dz = (Zmax - Zmin) / nSlicesZ;
178 Double_t dy = (Ymax - Ymin) / nSlicesY;
179
180 Int_t counter = 0;
181
182 for (Int_t i = 0; i < dimGrid; i++) {
183 if (fDebug)
184 if (i % 10 == 0)
185 cout << "Processing node " << i << endl;
186 bField[i] = new TGraph*[dimGrid];
187 Double_t ySlice = Ymin + i * dy;
188
189 for (Int_t j = 0; j < dimGrid; j++) {
190 Float_t zSlice = Zmin + j * dz;
191 bField[i][j] = new TGraph();
192
193 for (Int_t k = 0; k < fMap->GetNx(); k++) {
194 Double_t x = fMap->GetXmin() + k * fMap->GetXstep();
195 bField[i][j]->SetPoint(k, x,
196 (comp == "x") ? fMap->GetBx(x, ySlice, zSlice)
197 : (comp == "y") ? fMap->GetBy(x, ySlice, zSlice)
198 : (comp == "z") ? fMap->GetBz(x, ySlice, zSlice)
199 : 0.0);
200 }
201
202 vector<Double_t> fit_par;
203
204 TF1* f = nullptr;
205 if (iPart == 0) {
206 f = new TF1("f", "gaus", fConditions->tailRangeLeft.first, fConditions->tailRangeLeft.second);
207 f->SetParameter(2, 50.);
208 fit_par = bField[i][j]->Fit(f, "SQ", "", f->GetXmin(), f->GetXmax())->Parameters();
209 } else {
210 f = new TF1("f", "gaus", fConditions->tailRangeRight.first, fConditions->tailRangeRight.second);
211 f->SetParameter(2, 100.);
212 fit_par = bField[i][j]->Fit(f, "SQ", "", f->GetXmin(), f->GetXmax())->Parameters();
213 }
214
215 fMonitor->C[iPart][iComp]->Fill(fit_par[0]);
216 fMonitor->M[iPart][iComp]->Fill(fit_par[1]);
217 fMonitor->S[iPart][iComp]->Fill(fit_par[2]);
218
219 constant[iPart][iComp]->SetPoint(counter, ySlice, zSlice, fit_par[0]);
220 mean[iPart][iComp]->SetPoint(counter, ySlice, zSlice, fit_par[1]);
221 sigma[iPart][iComp]->SetPoint(counter, ySlice, zSlice, fit_par[2]);
222
223 if (f)
224 delete f;
225
226 counter++;
227 }
228 }
229
230 for (Int_t i = 0; i < dimNhist; i++) {
231 fit_func[i] = new TF1*[dimNhist];
232 Double_t y = WinYmin + i * winYstep;
233
234 for (Int_t j = 0; j < dimNhist; j++) {
235 Double_t z = WinZmin + j * winZstep;
236 fit_func[i][j] = new TF1("func", "[0]*exp(-0.5*((x-[1])/[2])**2)",
237 (iPart == 1) ? fConditions->tailRangeRight.second : fConditions->min,
238 (iPart == 1) ? fConditions->max : fConditions->tailRangeLeft.first);
239 fit_func[i][j]->SetParameter(0, constant[iPart][iComp]->Interpolate(y, z));
240 fit_func[i][j]->SetParameter(1, mean[iPart][iComp]->Interpolate(y, z));
241 fit_func[i][j]->SetParameter(2, sigma[iPart][iComp]->Interpolate(y, z));
242 }
243 }
244
245 DrawProfileExtrapX(dimNhist, bField, fit_func, "B" + comp + ".png", comp, Ymin, Zmin, dy, dz);
246
247 for (Int_t i = 0; i < dimNhist; i++) {
248 for (Int_t j = 0; j < dimNhist; j++)
249 delete fit_func[i][j];
250 delete[] fit_func[i];
251 }
252 delete[] fit_func;
253
254 for (Int_t i = 0; i < dimGrid; i++) {
255 for (Int_t j = 0; j < dimGrid; j++)
256 delete bField[i][j];
257 delete[] bField[i];
258 }
259 delete[] bField;
260}
261
262void BmnSP41FieldMapCreator::DoExtrapolationAlongZ(Int_t iPart, Int_t iComp)
263{
264 TGraph*** bField = new TGraph**[dimGrid];
265 TF1*** fit_func = new TF1**[dimNhist];
266 TString comp = (iComp == 0) ? "x" : (iComp == 1) ? "y" : (iComp == 2) ? "z" : "";
267
268 Double_t Xmax = fMap->GetXmax();
269 Double_t Xmin = fMap->GetXmin();
270
271 Double_t Ymax = fMap->GetYmax();
272 Double_t Ymin = fMap->GetYmin();
273
274 // Define number of slices on X(Y)-dimension of the grid used
275 Int_t nSlicesX = dimGrid - 1;
276 Int_t nSlicesY = dimGrid - 1;
277
278 // Define step along X(Y, Z)-axises of the used grid
279 Double_t dx = (Xmax - Xmin) / nSlicesX;
280 Double_t dy = (Ymax - Ymin) / nSlicesY;
281
282 Int_t counter = 0;
283
284 for (Int_t i = 0; i < dimGrid; i++) {
285 if (fDebug)
286 if (i % 10 == 0)
287 cout << "Processing node " << i << endl;
288 bField[i] = new TGraph*[dimGrid];
289 Double_t xSlice = Xmin + i * dx;
290
291 for (Int_t j = 0; j < dimGrid; j++) {
292 Float_t ySlice = Ymin + j * dy;
293 bField[i][j] = new TGraph();
294
295 for (Int_t k = 0; k < fMap->GetNz(); k++) {
296 Double_t z = fMap->GetZmin() + k * fMap->GetZstep();
297 bField[i][j]->SetPoint(k, z,
298 (comp == "x") ? fMap->GetBx(xSlice, ySlice, z)
299 : (comp == "y") ? fMap->GetBy(xSlice, ySlice, z)
300 : (comp == "z") ? fMap->GetBz(xSlice, ySlice, z)
301 : 0.0);
302 }
303
304 vector<Double_t> fit_par;
305
306 TF1* f = nullptr;
307 if (iPart == 0) {
308 f = new TF1("f", "gaus", fConditions->tailRangeLeft.first, fConditions->tailRangeLeft.second);
309 f->SetParameter(2, 100.);
310 fit_par = bField[i][j]->Fit(f, "SQ", "G", f->GetXmin(), f->GetXmax())->Parameters();
311 } else {
312 f = new TF1("f", "gaus", fConditions->tailRangeRight.first, fConditions->tailRangeRight.second);
313 f->SetParameter(2, 120.);
314 fit_par = bField[i][j]->Fit(f, "SQ", "G", f->GetXmin(), f->GetXmax())->Parameters();
315 }
316
317 fMonitor->C[iPart][iComp]->Fill(fit_par[0]);
318 fMonitor->M[iPart][iComp]->Fill(fit_par[1]);
319 fMonitor->S[iPart][iComp]->Fill(fit_par[2]);
320
321 constant[iPart][iComp]->SetPoint(counter, xSlice, ySlice, fit_par[0]);
322 mean[iPart][iComp]->SetPoint(counter, xSlice, ySlice, fit_par[1]);
323 sigma[iPart][iComp]->SetPoint(counter, xSlice, ySlice, fit_par[2]);
324
325 counter++;
326
327 if (f)
328 delete f;
329 }
330 }
331
332 for (Int_t i = 0; i < dimNhist; i++) {
333 fit_func[i] = new TF1*[dimNhist];
334 Double_t x = WinXmin + i * winXstep;
335
336 for (Int_t j = 0; j < dimNhist; j++) {
337 Double_t y = WinYmin + j * winYstep;
338 fit_func[i][j] = new TF1("func", "[0]*exp(-0.5*((x-[1])/[2])**2)",
339 (iPart == 1) ? fConditions->tailRangeRight.second : fConditions->min,
340 (iPart == 1) ? fConditions->max : fConditions->tailRangeLeft.first);
341 fit_func[i][j]->SetParameter(0, constant[iPart][iComp]->Interpolate(x, y));
342 fit_func[i][j]->SetParameter(1, mean[iPart][iComp]->Interpolate(x, y));
343 fit_func[i][j]->SetParameter(2, sigma[iPart][iComp]->Interpolate(x, y));
344 }
345 }
346
347 DrawProfileExtrapZ(dimNhist, bField, fit_func, "B" + comp + ".pdf", comp, Xmin, Ymin, dx, dy);
348
349 for (Int_t i = 0; i < dimNhist; i++) {
350 for (Int_t j = 0; j < dimNhist; j++)
351 delete fit_func[i][j];
352 delete[] fit_func[i];
353 }
354 delete[] fit_func;
355
356 for (Int_t i = 0; i < dimGrid; i++) {
357 for (Int_t j = 0; j < dimGrid; j++)
358 delete bField[i][j];
359 delete[] bField[i];
360 }
361 delete[] bField;
362}
363
364void BmnSP41FieldMapCreator::DrawProfileExtrapZ(Int_t dim,
365 TGraph*** gr,
366 TF1*** f,
367 TString name,
368 TString comp,
369 Double_t Xmin,
370 Double_t Ymin,
371 Double_t dx,
372 Double_t dy)
373{
374 TCanvas* c = new TCanvas("1", "1", 1200, 800);
375 c->Divide(dim, dim);
376 Int_t counter = 1;
377 for (Int_t i = 0; i < dim; i++) {
378 Double_t X = WinXmin + i * winXstep;
379 Double_t sliceX = (X - Xmin) / dx;
380 Int_t shiftI = Int_t(sliceX);
381 for (Int_t j = 0; j < dim; j++) {
382 Double_t Y = WinYmin + j * winYstep;
383 Double_t sliceY = (Y - Ymin) / dy;
384 Int_t shiftJ = Int_t(sliceY);
385 c->cd(counter);
386 gr[shiftI][shiftJ]->Draw();
387 gr[shiftI][shiftJ]->GetXaxis()->SetLimits(fConditions->min, fConditions->max);
388 gr[shiftI][shiftJ]->SetTitle("B" + comp + TString(Form(", x = %G [cm]", X))
389 + TString(Form(", y = %G [cm]", Y)));
390 gr[shiftI][shiftJ]->GetXaxis()->SetTitleSize(0.1);
391 gr[shiftI][shiftJ]->GetXaxis()->SetLabelSize(0.08);
392 gr[shiftI][shiftJ]->GetYaxis()->SetTitleSize(0.1);
393 gr[shiftI][shiftJ]->GetYaxis()->SetLabelSize(0.1);
394 gr[shiftI][shiftJ]->GetXaxis()->SetTitle("Z, cm");
395 gr[shiftI][shiftJ]->GetYaxis()->SetTitle("B" + comp + ", kG");
396 counter++;
397 f[i][j]->Draw("same");
398 }
399 }
400 c->SaveAs(name);
401 delete c;
402}
403
404void BmnSP41FieldMapCreator::DrawProfileExtrapX(Int_t dim,
405 TGraph*** gr,
406 TF1*** f,
407 TString name,
408 TString comp,
409 Double_t Ymin,
410 Double_t Zmin,
411 Double_t dy,
412 Double_t dz)
413{
414 gStyle->SetTitleFontSize(0.1);
415 TCanvas* c = new TCanvas("1", "1", 4000, 4000);
416 c->Divide(dim, dim);
417 Int_t counter = 1;
418 for (Int_t i = 0; i < dim; i++) {
419 Double_t Y = WinYmin + i * winYstep;
420 Double_t sliceY = (Y - Ymin) / dy;
421 Int_t shiftI = Int_t(sliceY);
422 for (Int_t j = 0; j < dim; j++) {
423 Double_t Z = WinZmin + j * winZstep;
424 Double_t sliceZ = (Z - Zmin) / dz;
425 Int_t shiftJ = Int_t(sliceZ);
426 c->cd(counter);
427 gr[shiftI][shiftJ]->Draw();
428 gr[shiftI][shiftJ]->GetXaxis()->SetLimits(fConditions->min, fConditions->max);
429 gr[shiftI][shiftJ]->SetLineWidth(2);
430 gr[shiftI][shiftJ]->GetXaxis()->SetNdivisions(5);
431 gr[shiftI][shiftJ]->SetTitle("B" + comp + TString(Form(", y = %G [cm]", Y))
432 + TString(Form(", z = %G [cm]", Z)));
433 gr[shiftI][shiftJ]->GetXaxis()->SetTitleSize(0.1);
434 gr[shiftI][shiftJ]->GetXaxis()->SetLabelSize(0.08);
435 gr[shiftI][shiftJ]->GetYaxis()->SetTitleSize(0.1);
436 gr[shiftI][shiftJ]->GetYaxis()->SetLabelSize(0.1);
437 gr[shiftI][shiftJ]->GetXaxis()->SetTitle("X, cm");
438 gr[shiftI][shiftJ]->GetYaxis()->SetTitle("B" + comp + ", kG");
439 counter++;
440 f[i][j]->Draw("same");
441 f[i][j]->SetLineWidth(2);
442 }
443 }
444 c->SaveAs(name);
445 delete c;
446}
447
448void BmnSP41FieldMapCreator::CreateExtrapolatedMap()
449{
450
451 Double_t Xmin = fMap->GetXmin();
452 Double_t Ymin = fMap->GetYmin();
453 Double_t Zmin = fMap->GetZmin();
454
455 Double_t Xmax = fMap->GetXmax();
456 Double_t Ymax = fMap->GetYmax();
457 Double_t Zmax = fMap->GetZmax();
458
459 Double_t dx = fMap->GetXstep();
460 Double_t dy = fMap->GetYstep();
461 Double_t dz = fMap->GetZstep();
462
463 Int_t Nx = fMap->GetNx();
464 Int_t Ny = fMap->GetNy();
465 Int_t Nz = fMap->GetNz();
466 // Int_t Nz = (Zmax - Zmin) / dz + 1;
467
468 // Print(fMap);
469
470 Double_t x = -1., y = -1., z = -1.;
471
472 ofstream mapFile(nameOutput);
473
474 mapFile.precision(6);
475 mapFile << showpoint;
476 mapFile << "nosym" << endl;
477
478 TString axis = fConditions->direction;
479
480 // Filling parameters of the extrapolated map ...
481 if (axis.Contains("X")) {
482 Xmin = fConditions->min;
483 Xmax = fConditions->max;
484 Nx = Int_t((Xmax - Xmin) / 2.) + 1;
485 mapFile << Xmin << " " << Xmax << " " << Nx << endl;
486 } else
487 mapFile << Xmin << " " << Xmax << " " << Nx << endl;
488
489 mapFile << Ymin << " " << Ymax << " " << Ny << endl;
490
491 if (axis.Contains("Z")) {
492 Zmin = fConditions->min;
493 Zmax = fConditions->max;
494 Nz = Int_t((Zmax - Zmin) / 2.) + 1;
495 mapFile << Zmin << " " << Zmax << " " << Nz << endl;
496 } else
497 mapFile << Zmin << " " << Zmax << " " << Nz << endl;
498
499 for (Int_t ix = 0; ix < Nx; ix++) {
500 x = Xmin + dx * ix;
501 for (Int_t iy = 0; iy < Ny; iy++) {
502 y = Ymin + dy * iy;
503 if (fDebug)
504 cout << x << " " << y << endl;
505 for (Int_t iz = 0; iz < Nz; iz++) {
506 // Int_t ind = ix * Ny * Nz + iy * Nz + iz;
507 z = Zmin + dz * iz;
508
509 Double_t bx = -1., by = -1., bz = -1.;
510
511 Bool_t isInMap = kTRUE;
512
513 if (axis.Contains("X")) {
514 if (x < fMap->GetXmin()) {
515 isInMap = kFALSE;
516 bx = FieldExtrapolate(x, y, z, constant[0][0], mean[0][0], sigma[0][0]);
517 by = FieldExtrapolate(x, y, z, constant[0][1], mean[0][1], sigma[0][1]);
518 bz = FieldExtrapolate(x, y, z, constant[0][2], mean[0][2], sigma[0][2]);
519 } else if (x > fMap->GetXmax() - 1.) {
520 isInMap = kFALSE;
521 bx = FieldExtrapolate(x, y, z, constant[1][0], mean[1][0], sigma[1][0]);
522 by = FieldExtrapolate(x, y, z, constant[1][1], mean[1][1], sigma[1][1]);
523 bz = FieldExtrapolate(x, y, z, constant[1][2], mean[1][2], sigma[1][2]);
524 }
525 } else if (axis.Contains("Z")) {
526 if (z < fMap->GetZmin()) {
527 isInMap = kFALSE;
528 bx = FieldExtrapolate(x, y, z, constant[0][0], mean[0][0], sigma[0][0]);
529 by = FieldExtrapolate(x, y, z, constant[0][1], mean[0][1], sigma[0][1]);
530 bz = FieldExtrapolate(x, y, z, constant[0][2], mean[0][2], sigma[0][2]);
531 } else if (z > fMap->GetZmax() - 1.) {
532 isInMap = kFALSE;
533 bx = FieldExtrapolate(x, y, z, constant[1][0], mean[1][0], sigma[1][0]);
534 by = FieldExtrapolate(x, y, z, constant[1][1], mean[1][1], sigma[1][1]);
535 bz = FieldExtrapolate(x, y, z, constant[1][2], mean[1][2], sigma[1][2]);
536 }
537 }
538
539 if (isInMap) {
540 bx = fMap->GetBx(x, y, z);
541 by = fMap->GetBy(x, y, z);
542 bz = fMap->GetBz(x, y, z);
543 }
544
545 // Checking obtained extrapolated values ...
546 mapFile << checkedFieldValue(bx, by, bz).at(0) / factor << " "
547 << checkedFieldValue(bx, by, bz).at(1) / factor << " "
548 << checkedFieldValue(bx, by, bz).at(2) / factor << endl;
549 } // z-Loop
550 } // y-Loop
551 } // x-Loop
552}
553
554void BmnSP41FieldMapCreator::SmoothMap(TClonesArray* sValues)
555{
556
557 vector<TLorentzVector> vSmoo;
558
559 for (Int_t iSmooth = 0; iSmooth < sValues->GetEntriesFast(); iSmooth++) {
560 smoothedValues* sValue = (smoothedValues*)sValues->UncheckedAt(iSmooth);
561
562 TLorentzVector tmp(sValue->X, sValue->Y, sValue->Z, sValue->fieldValue);
563 vSmoo.push_back(tmp);
564 }
565
566 Double_t Xmin = fMap->GetXmin();
567 Double_t Ymin = fMap->GetYmin();
568 Double_t Zmin = fMap->GetZmin();
569 Double_t Zmax = fMap->GetZmax();
570
571 Double_t dx = fMap->GetXstep();
572 Double_t dy = fMap->GetYstep();
573 Double_t dz = fMap->GetZstep();
574
575 Int_t Nx = fMap->GetNx();
576 Int_t Ny = fMap->GetNy();
577 Int_t Nz = fMap->GetNz();
578
579 Print(fMap);
580
581 Double_t x, y, z, bx, by, bz;
582
583 UInt_t index = 0;
584 ofstream mapFile(nameOutput);
585
586 mapFile.precision(6);
587 mapFile << showpoint;
588 mapFile << "nosym" << endl;
589 mapFile << Xmin << " " << -Xmin << " " << Nx << endl;
590 mapFile << Ymin << " " << -Ymin << " " << Ny << endl;
591 mapFile << Zmin << " " << Zmax << " " << Nz << endl;
592
593 for (Int_t ix = 0; ix < Nx; ix++) {
594 x = Xmin + dx * ix;
595
596 for (Int_t iy = 0; iy < Ny; iy++) {
597 y = Ymin + dy * iy;
598
599 for (Int_t iz = 0; iz < Nz; iz++) {
600 index = ix * Ny * Nz + iy * Nz + iz;
601 z = Zmin + dz * iz;
602
603 by = (fMap->GetBy(x, y, z) / factor);
604 for (auto it = vSmoo.begin(); it != vSmoo.end(); it++) {
605
606 TLorentzVector v = *it;
607
608 if (v.X() != x || v.Y() != y || v.Z() != z)
609 continue;
610
611 by = v.T() / factor;
612
613 vSmoo.erase(it--);
614 break;
615 }
616
617 bx = (Abs(fMap->GetBx(x, y, z)) / factor > 10e-8) ? (fMap->GetBx(x, y, z) / factor) : 0.0;
618 bz = (Abs(fMap->GetBz(x, y, z)) / factor > 10e-8) ? (fMap->GetBz(x, y, z) / factor) : 0.0;
619 // mapFile << x << " " << y << " " << z << " " << bx << " " << by << " " << bz << endl;
620 mapFile << bx << " " << by << " " << bz << endl;
621
622 if (index % 1000 == 0)
623 cout << "nNodes processed# " << index << endl;
624
625 } // z-Loop
626 } // y-Loop
627 } // x-Loop
628}
629
630Double_t BmnSP41FieldMapCreator::FieldExtrapolate(Double_t x,
631 Double_t y,
632 Double_t z,
633 TGraph2D* c,
634 TGraph2D* m,
635 TGraph2D* s)
636{
637 vector<Double_t> varToBePassed;
638
639 if (fConditions->direction.Contains("X")) {
640 varToBePassed.push_back(y);
641 varToBePassed.push_back(z);
642 varToBePassed.push_back(x);
643 } else if (fConditions->direction.Contains("Z")) {
644 varToBePassed.push_back(x);
645 varToBePassed.push_back(y);
646 varToBePassed.push_back(z);
647 }
648
649 Double_t Mean = m->Interpolate(varToBePassed[0], varToBePassed[1]);
650 Double_t Sigma = s->Interpolate(varToBePassed[0], varToBePassed[1]);
651
652 return c->Interpolate(varToBePassed[0], varToBePassed[1])
653 * TMath::Exp(-0.5 * (varToBePassed[2] - Mean) * (varToBePassed[2] - Mean) / Sigma / Sigma);
654}
655
656void BmnSP41FieldMapCreator::Print(BmnFieldMap* magField)
657{
658 cout << "Nx = " << magField->GetNx();
659 cout << "; Ny = " << magField->GetNy();
660 cout << "; Nz = " << magField->GetNz();
661
662 Double_t shift = 1e-6;
663
664 cout << "; Xmin = " << magField->GetXmin() - shift;
665 cout << "; dx = " << magField->GetXstep();
666 cout << "; Ymin = " << magField->GetYmin() - shift;
667 cout << "; dy = " << magField->GetYstep();
668 cout << "; Zmin = " << magField->GetZmin() - shift;
669 cout << "; dz = " << magField->GetZstep() << endl;
670}
v
__m128 v
Definition P4_F32vec4.h:1
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
Bool_t
bool Bool_t
Int_t
int Int_t
Float_t
float Float_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Form
char * Form(const char *fmt,...)
gStyle
R__EXTERN TStyle * gStyle
BmnFieldMap::GetNz
Int_t GetNz() const
Definition BmnFieldMap.h:82
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::GetZmin
Double_t GetZmin() const
Definition BmnFieldMap.h:70
BmnFieldMap::GetYmin
Double_t GetYmin() const
Definition BmnFieldMap.h:69
BmnFieldMap::GetNy
Int_t GetNy() const
Definition BmnFieldMap.h:81
BmnFieldMap::GetNx
Int_t GetNx() const
Definition BmnFieldMap.h:80
BmnFieldMap::GetXmin
Double_t GetXmin() const
Accessors to field parameters in local coordinate system.
Definition BmnFieldMap.h:68
BmnNewFieldMap
Field Map does not have any symmetries w.r.t coordinate axises.
Definition BmnNewFieldMap.h:16
BmnSP41FieldMapCreator::FieldExtrapolate
Double_t FieldExtrapolate(Double_t, Double_t, Double_t, TGraph2D *, TGraph2D *, TGraph2D *)
Definition BmnSP41FieldMapCreator.cxx:630
BmnSP41FieldMapCreator::DrawProfileExtrapZ
void DrawProfileExtrapZ(Int_t, TGraph ***, TF1 ***, TString, TString, Double_t, Double_t, Double_t, Double_t)
Definition BmnSP41FieldMapCreator.cxx:364
BmnSP41FieldMapCreator::checkedFieldValue
std::vector< Double_t > checkedFieldValue(Double_t bx, Double_t by, Double_t bz)
Definition BmnSP41FieldMapCreator.h:145
BmnSP41FieldMapCreator::DoExtrapolationAlongZ
void DoExtrapolationAlongZ(Int_t, Int_t)
Definition BmnSP41FieldMapCreator.cxx:262
BmnSP41FieldMapCreator::DrawProfileExtrapX
void DrawProfileExtrapX(Int_t, TGraph ***, TF1 ***, TString, TString, Double_t, Double_t, Double_t, Double_t)
Definition BmnSP41FieldMapCreator.cxx:404
BmnSP41FieldMapCreator::DoExtrapolationAlongX
void DoExtrapolationAlongX(Int_t, Int_t)
Definition BmnSP41FieldMapCreator.cxx:160
BmnSP41FieldMapCreator::fConditions
extrapolationConditions * fConditions
Definition BmnSP41FieldMapCreator.h:118
ROOT::Fit::FitResult::Parameters
const std::vector< double > & Parameters() const
TAttAxis::SetLabelSize
virtual void SetLabelSize(Float_t size=0.04)
TAttAxis::SetTitleSize
virtual void SetTitleSize(Float_t size=0.04)
TAttAxis::SetNdivisions
virtual void SetNdivisions(Int_t n1, Int_t n2, Int_t n3, Bool_t optim=kTRUE)
TAttLine::SetLineWidth
virtual void SetLineWidth(Width_t lwidth)
TAxis::SetLimits
virtual void SetLimits(Double_t xmin, Double_t xmax)
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
TF1::SetParameter
virtual void SetParameter(const TString &name, Double_t value)
TGraph2D::Interpolate
Double_t Interpolate(Double_t x, Double_t y)
TGraph::SetPoint
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
TGraph::Fit
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Axis_t xmin=0, Axis_t xmax=0)
TGraph::Draw
void Draw(Option_t *chopt="") override
TGraph::GetXaxis
TAxis * GetXaxis() const
TGraph::GetYaxis
TAxis * GetYaxis() const
TGraph::SetTitle
void SetTitle(const char *title="") override
TNamed::SetTitle
virtual void SetTitle(const char *title="")
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::UncheckedAt
TObject * UncheckedAt(Int_t i) const
TObject::TClonesArray
friend friend class TClonesArray
TPad::Divide
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
TPad::SaveAs
void SaveAs(const char *filename="", Option_t *option="") const override
TString::Data
const char * Data() const
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
dim
const UInt_t dim
Definition MpdVHLLEGenerator.h:22
TMath::Exp
Double_t Exp(Double_t x)
TMath::Sigma
constexpr Double_t Sigma()
TMath::Abs
Double_t Abs(Double_t d)
TMath::Mean
Double_t Mean(Iterator first, Iterator last)
std
Definition json.hpp:4446