29 #include "FairRunAna.h"
30 #include "FairField.h"
39 : fNStations(0), fStations(0)
48 float xMax[48] = { 66, 66, 66, 66, 66, 66, 66, 66,
49 66, 66, 66, 66, 66, 66, 66, 66,
50 88, 88, 88, 88, 88, 88, 88, 88,
51 104,104,104,104,104,104,104,104,
52 160,160,160,160,160,160,160,160,
53 160,160,160,160,160,160,160,160 };
55 float yMax[48] = { 32, 32, 32, 32, 32, 32, 32, 32,
56 32, 32, 32, 32, 32, 32, 32, 32,
57 37, 37, 37, 37, 37, 37, 37, 37,
58 41, 41, 41, 41, 41, 41, 41, 41,
59 60, 60, 60, 60, 60, 60, 60, 60,
60 60, 60, 60, 60, 60, 60, 60, 60 };
66 if( ist == 32 || ist == 16)
68 vector<L1FieldSlice> &virtualSlices =
fStations[ist].fieldVirtualSlice;
69 virtualSlices.resize(10);
71 float xMaxVirtual = xMax[ist];
72 float yMaxVirtual = yMax[ist];
73 float dx = (xMax[ist] - xMax[ist-1])/11.
f;
74 float dy = (yMax[ist] - yMax[ist-1])/11.
f;
76 for(
unsigned int iV = 0; iV < virtualSlices.size(); iV++ )
99 if( dx > xMax/N/2 ) dx = xMax/N/4.;
100 if( dy > yMax/N/2 ) dy = yMax/N/4.;
103 TVectorD b0(N),
b1(N),
b2(N);
104 for(
int i=0;
i<N;
i++){
105 for(
int j=0; j<N; j++) A(
i,j) = 0.;
109 if(FairRunAna::Instance()->GetField())
111 for(
double x=-xMax;
x<=xMax;
x+=
dx )
112 for(
double y=-yMax;
y<=yMax;
y+=
dy )
119 FairRunAna::Instance()->GetField()->Field(p, B);
122 for(
int i=1;
i<=M;
i++){
125 for(
int j=0; j<
i; j++ )
m(l+j) =
x*
m(k+j);
130 for(
int i=0;
i<N;
i++){
131 for(
int j=0; j<N;j++) A(
i,j)+=w*
m(
i)*
m(j);
139 TVectorD c0 = A*b0,
c1 = A*
b1,
c2 = A*
b2;
140 for(
int i=0;
i<N;
i++){
141 fieldSlice.
cx[
i] = c0(
i);
150 for(
int iV=0; iV<int_v::Size; iV++)
152 if(!(mask[iV]))
continue;
172 cout <<
"ERROR: Settings.data format is wrong! Station: " <<
i-1 << endl;
182 p.
fStations[
i].materialInfo.thick = xTimesRho/1.39;
316 TDirectory *curr = gDirectory;
317 TFile *currentFile = gFile;
318 TFile* fout =
new TFile(
"FieldApprox.root",
"RECREATE");
321 float xMax[48] = { 66, 66, 66, 66, 66, 66, 66, 66,
322 66, 66, 66, 66, 66, 66, 66, 66,
323 88, 88, 88, 88, 88, 88, 88, 88,
324 104,104,104,104,104,104,104,104,
325 160,160,160,160,160,160,160,160,
326 160,160,160,160,160,160,160,160 };
328 float yMax[48] = { 32, 32, 32, 32, 32, 32, 32, 32,
329 32, 32, 32, 32, 32, 32, 32, 32,
330 37, 37, 37, 37, 37, 37, 37, 37,
331 41, 41, 41, 41, 41, 41, 41, 41,
332 60, 60, 60, 60, 60, 60, 60, 60,
333 60, 60, 60, 60, 60, 60, 60, 60 };
338 double Xmax=xMax[ist], Ymax=yMax[ist];
348 float ddx = 2*Xmax/NbinsX;
349 float ddy = 2*Ymax/NbinsY;
351 TH2F *stB =
new TH2F(Form(
"station %i, dB", ist+1) ,Form(
"station %i, dB, z = %0.f cm", ist+1,z) , static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
352 TH2F *stBx =
new TH2F(Form(
"station %i, dBx", ist+1),Form(
"station %i, dBx, z = %0.f cm", ist+1,z), static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
353 TH2F *stBy =
new TH2F(Form(
"station %i, dBy", ist+1),Form(
"station %i, dBy, z = %0.f cm", ist+1,z), static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
354 TH2F *stBz =
new TH2F(Form(
"station %i, dBz", ist+1),Form(
"station %i, dBz, z = %0.f cm", ist+1,z), static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
361 for(
int ii = 1; ii <=NbinsX+1; ii++)
364 x = -Xmax+(ii-1)*ddx;
365 for(
int jj = 1;
jj <=NbinsY+1;
jj++)
367 y = -Ymax+(
jj-1)*ddy;
374 r[2] =
z; r[0] =
x; r[1] =
y;
375 FairRunAna::Instance()->GetField()->Field(r, B);
376 bbb =
sqrt(B[0]*B[0]+B[1]*B[1]+B[2]*B[2]);
379 bbb_L1 =
sqrt(B_L1.x[0]*B_L1.x[0] + B_L1.y[0]*B_L1.y[0] + B_L1.z[0]*B_L1.z[0]);
381 stB -> SetBinContent(ii,
jj,
fabs(bbb)>1.e-2 ? (bbb - bbb_L1)/bbb*100 : 0);
382 stBx -> SetBinContent(ii,
jj,
fabs(B[0])>1.e-2 ? (B[0] - B_L1.x[0])/B[0]*100 : 0);
383 stBy -> SetBinContent(ii,
jj,
fabs(B[1])>1.e-2 ? (B[1] - B_L1.y[0])/B[1]*100 : 0);
384 stBz -> SetBinContent(ii,
jj,
fabs(B[2])>1.e-2 ? (B[2] - B_L1.z[0])/B[2]*100 : 0);
400 stB ->GetXaxis()->SetTitle(
"X, cm");
401 stB ->GetYaxis()->SetTitle(
"Y, cm");
402 stB ->GetXaxis()->SetTitleOffset(1);
403 stB ->GetYaxis()->SetTitleOffset(1);
404 stB ->GetZaxis()->SetTitle(
"B_map - B_L1, kGauss");
405 stB ->GetZaxis()->SetTitleOffset(1.3);
407 stBx ->GetXaxis()->SetTitle(
"X, cm");
408 stBx ->GetYaxis()->SetTitle(
"Y, cm");
409 stBx ->GetXaxis()->SetTitleOffset(1);
410 stBx ->GetYaxis()->SetTitleOffset(1);
411 stBx ->GetZaxis()->SetTitle(
"Bx_map - Bx_L1, kGauss");
412 stBx ->GetZaxis()->SetTitleOffset(1.3);
414 stBy ->GetXaxis()->SetTitle(
"X, cm");
415 stBy ->GetYaxis()->SetTitle(
"Y, cm");
416 stBy ->GetXaxis()->SetTitleOffset(1);
417 stBy ->GetYaxis()->SetTitleOffset(1);
418 stBy ->GetZaxis()->SetTitle(
"By_map - By_L1, kGauss");
419 stBy ->GetZaxis()->SetTitleOffset(1.3);
421 stBz ->GetXaxis()->SetTitle(
"X, cm");
422 stBz ->GetYaxis()->SetTitle(
"Y, cm");
423 stBz ->GetXaxis()->SetTitleOffset(1);
424 stBz ->GetYaxis()->SetTitleOffset(1);
425 stBz ->GetZaxis()->SetTitle(
"Bz_map - Bz_L1, kGauss");
426 stBz ->GetZaxis()->SetTitleOffset(1.3);
friend F32vec4 cos(const F32vec4 &a)
const FTSCAStation & Station(short i) const
void GetStripInfo(FTSCAStripInfoVector &stripInfo, const int_v iStation, const float_m &mask) const
friend F32vec4 sqrt(const F32vec4 &a)
friend F32vec4 sin(const F32vec4 &a)
friend F32vec4 log(const F32vec4 &a)
void GetFieldValue(const float_v &x, const float_v &y, CAFieldValue &B, const float_m &mask=float_m(true)) const
static void BinaryStoreWrite(const T *mem, int count, FILE *f)
CAFieldValue fVtxFieldValue[2]
friend F32vec4 fabs(const F32vec4 &a)
static void BinaryStoreRead(T *mem, int count, FILE *f)
void RestoreFromFile(FILE *f)
void CalculateFieldSlice(L1FieldSlice &fieldSlice, const float xMax, const float yMax, const float z)
void CheckFieldApproximation()
TBuffer & operator>>(TBuffer &buf, PndAnaPidSelector *&obj)
void StoreToFile(FILE *f) const
TMatrixT< double > TMatrixD