]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCEfield.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
remove printouts (Sandro B.)
[u/mrichter/AliRoot.git] / TPC / AliTPCEfield.cxx
CommitLineData
b1f0a2a5 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17 Calculation of the Electric field:
18 Sollution of laplace equation in cartezian system, with boundary condition.
19 Se details:
20 http://web.mit.edu/6.013_book/www/chapter5/5.10.html
21
22
23*/
24
25/* $Id: AliTPCEfield.cxx 41275 2010-05-16 22:23:06Z marian $ */
26
27#include "TTreeStream.h"
28#include "TMath.h"
29#include "TLinearFitter.h"
30#include "TRandom.h"
31#include "AliTPCEfield.h"
32
33ClassImp(AliTPCEfield)
34
35
36AliTPCEfield* AliTPCEfield::fgInstance=0;
37
38AliTPCEfield::AliTPCEfield():
39 TNamed(),
40 fScale(0),
41 fMaxFreq(0),
42 fIs2D(kTRUE),
43 fWorkspace(0) // file with trees, pictures ...
44{
45 //
46 for (Int_t i=0; i<3; i++){
47 fMin[i]=0; fMax[i]=0;
48 }
49 fgInstance =this;
50}
51
52AliTPCEfield::AliTPCEfield(const char* name, Int_t maxFreq, Bool_t is2D, Bool_t useLinear):
53 TNamed(name,name),
54 fScale(0),
55 fMaxFreq(maxFreq),
56 fIs2D(is2D),
57 fUseLinear(useLinear),
58 fWorkspace(0) // file with trees, pictures ...
59{
60 //
61 for (Int_t i=0; i<3; i++){
62 fMin[i]=0; fMax[i]=0;
63 }
64 fWorkspace=new TTreeSRedirector(Form("%s.root",name));
65 MakeFitFunctions(maxFreq);
66 fgInstance =this;
67}
68
69void AliTPCEfield::MakeFitFunctions(Int_t maxFreq){
70 //
71 // fit functions = f(x,y,z) = fx(x)*fy(y)*fz(z)
72 // function can be of following types:
73 // 0 - constant
74 // 1 - linear
75 // 2 - hypx
76 // 3 - hypy
77 // 4 - hypz
78 //
79 Int_t nfunctions=0;
80 if (fIs2D) nfunctions = 1+(maxFreq)*8;
81 if (!fIs2D) nfunctions = 1+(maxFreq)*8;
82 if (fUseLinear) nfunctions+=2;
83 fFitFunctions = new TMatrixD(nfunctions,4);
84 fFitParam = new TVectorD(nfunctions);
85 fFitCovar = new TMatrixD(nfunctions,nfunctions);
86 TMatrixD &fitF = *fFitFunctions;
87 // constant function
88 Int_t counter=1;
89 fitF(0,0)=0;
90 //
91 // linear functions in one cordinate - constan in others
92 if (fUseLinear)
93 for (Int_t ifx=0; ifx<=1; ifx++)
94 for (Int_t ify=0; ify<=1; ify++)
95 for (Int_t ifz=0; ifz<=1; ifz++){
96 if (fIs2D && ifz>0) continue;
97 if ((ifx+ify+ifz)==0) continue;
98 if ((ifx+ify+ifz)>1) continue;
99 fitF(counter,0)= 1; // function type
100 fitF(counter,1)= ifx; // type of x function
101 fitF(counter,2)= ify; // type of y function
102 fitF(counter,3)= ifz; // type of z function
103 counter++;
104 }
105 //
106 if (fIs2D){
107 for (Int_t ihyp=0; ihyp<2; ihyp++)
108 for (Int_t ifx=-maxFreq; ifx<=maxFreq; ifx++)
109 for (Int_t ify=-maxFreq; ify<=maxFreq; ify++){
110 if (TMath::Abs(ify)!=TMath::Abs(ifx)) continue;
111 if (ifx==0) continue;
112 if (ify==0) continue;
113 fitF(counter,0)= 2+ihyp; // function type
114 fitF(counter,1)= ifx; // type of x function - + sinus - cosin
115 fitF(counter,2)= ify; // type of y function
116 fitF(counter,3)= 0; // type of y function
117 counter++;
118 }
119 }
120
121}
122
123
124
125AliTPCEfield::~AliTPCEfield() {
126 //
127 // Destructor
128 //
129 if (fWorkspace) delete fWorkspace;
130}
131
132void AliTPCEfield::SetRange(Double_t x0, Double_t x1, Double_t y0, Double_t y1, Double_t z0,Double_t z1){
133 //
134 // Set the ranges - coordinates are rescaled in order to use proper
135 // cos,sin expansion in scaled space
136 //
137 fMin[0]=x0; fMax[0]=x1;
138 fMin[1]=y0; fMax[1]=y1;
139 fMin[2]=z0; fMax[2]=z1;
140 if (fIs2D) fScale=0.5*(TMath::Abs(x1-x0)+TMath::Abs(y1-y0));
141 if (!fIs2D) fScale=0.5*(TMath::Abs(x1-x0)+TMath::Abs(y1-y0)+TMath::Abs(z1-z0));
142}
143
144
145void AliTPCEfield::AddBoundaryLine(Double_t x0,Double_t y0,Double_t z0, Double_t v0, Double_t x1, Double_t y1, Double_t z1,Double_t v1, Int_t id, Int_t npoints){
146 //
147 // Add a e field boundary line
148 // From point (x0,y0) to point (x1,y1)
149 // Linear decrease of potential is assumed
150 // Boundary can be identified using boundary ID
151 // The line is written into tree Boundary
152 //
153 Double_t deltaX = (x1-x0);
154 Double_t deltaY = (y1-y0);
155 Double_t deltaZ = (z1-z0);
156 Double_t deltaV = (v1-v0);
157 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
158 Double_t bpoint=gRandom->Rndm();
159 Double_t x = x0+deltaX*bpoint;
160 Double_t y = y0+deltaY*bpoint;
161 Double_t z = z0+deltaZ*bpoint;
162 Double_t v = v0+deltaV*bpoint;
163 (*fWorkspace)<<"Boundary"<<
164 "x="<<x<< // x coordinate
165 "y="<<y<< // y coordinate
166 "z="<<z<< // z coordinate
167 "v="<<v<< // potential
168 "id="<<id<< // boundary ID
169 "\n";
170 }
171}
172
173TTree * AliTPCEfield::GetTree(const char * tname){
174 //
175 //
176 //
177 return ((*fWorkspace)<<tname).GetTree();
178}
179
180Double_t AliTPCEfield::Field(Int_t ftype, Double_t ifx, Double_t ify, Double_t ifz, Double_t x, Double_t y, Double_t z){
181 //
182 // Field component in
183 // f frequency
184 Double_t fx=1,fy=1,fz=1;
185 const Double_t kEps=0.01;
186 //
187 if (ftype==0) return 1;
188 if (ftype==1) {
189 if (TMath::Nint(ifx)==1) return x;
190 if (TMath::Nint(ify)==1) return y;
191 if (TMath::Nint(ifz)==1) return z;
192 }
193 Double_t pi = TMath::Pi();
194 if (ifx>kEps) fx = (ftype==2) ? SinHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Sin(ifx*pi*x);
195 if (ifx<-kEps) fx = (ftype==2) ? CosHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Cos(ifx*pi*x);
196 //
197 if (ify>kEps) fy = (ftype==3) ? SinHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Sin(ify*pi*y);
198 if (ify<-kEps) fy = (ftype==3) ? CosHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Cos(ify*pi*y);
199 //
200 if (ifz>kEps) fz = (ftype==4) ? SinHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Sin(ifz*pi*z);
201 if (ifz<-kEps) fz = (ftype==4) ? CosHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Cos(ifz*pi*z);
202 (*fWorkspace)<<"eval"<<
203 "x="<<x<<
204 "y="<<y<<
205 "z="<<z<<
206 "ifx="<<ifx<<
207 "ify="<<ify<<
208 "ifz="<<ifz<<
209 "fx="<<fx<<
210 "fy="<<fy<<
211 "fz="<<fz<<
212 "ftype="<<ftype<<
213 "\n";
214 return fx*fy*fz;
215}
216
217
218Double_t AliTPCEfield::FieldDn(Int_t ftype, Double_t ifx, Double_t ify, Double_t ifz, Int_t dn, Double_t x, Double_t y, Double_t z){
219 //
220 //
221 //
222 //
223 // Field component in
224 // f frequency
225 Double_t fx=1,fy=1,fz=1;
226 const Double_t kEps=0.01;
227 //
228 if (ftype==0) return 0.;
229 if (ftype==1) {
230 Double_t value=0;
231 if (TMath::Nint(ifx)==1 &&dn==0) value=1.;
232 if (TMath::Nint(ify)==1 &&dn==1) value=1.;
233 if (TMath::Nint(ifz)==1 &&dn==2) value=1.;
234 return value;
235 }
236 Double_t pi = TMath::Pi();
237 if (ifx>kEps) fx = (ftype==2) ? SinHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Sin(ifx*pi*x);
238 if (ifx<-kEps) fx = (ftype==2) ? CosHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Cos(ifx*pi*x);
239 //
240 if (ify>kEps) fy = (ftype==3) ? SinHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Sin(ify*pi*y);
241 if (ify<-kEps) fy = (ftype==3) ? CosHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Cos(ify*pi*y);
242 //
243 if (ifz>kEps) fz = (ftype==4) ? SinHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Sin(ifz*pi*z);
244 if (ifz<-kEps) fz = (ftype==4) ? CosHNorm(ifz*pi*z,-TMath::Abs(ifz*pi)): TMath::Cos(ifz*pi*z);
245
246 if (dn==0){
247 if (ifx>kEps) fx = (ftype==2) ? CosHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): TMath::Cos(ifx*pi*x);
248 if (ifx<-kEps) fx = (ftype==2) ? SinHNorm(ifx*pi*x,TMath::Abs(ifx*pi)): -TMath::Sin(ifx*pi*x);
249 fx*=ifx*pi;
250 }
251 if (dn==1){
252 if (ify>kEps) fy = (ftype==3) ? CosHNorm(ify*pi*y,TMath::Abs(ify*pi)): TMath::Cos(ify*pi*y);
253 if (ify<-kEps) fy = (ftype==3) ? SinHNorm(ify*pi*y,TMath::Abs(ify*pi)): -TMath::Sin(ify*pi*y);
254 fy*=ify*pi;
255 }
256 if (dn==2){
257 if (ifz>kEps) fz = (ftype==4) ? CosHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): TMath::Cos(ifz*pi*z);
258 if (ifz<-kEps) fz = (ftype==4) ? SinHNorm(ifz*pi*z,TMath::Abs(ifz*pi)): -TMath::Sin(ifz*pi*z);
259 fz*=ifz*pi;
260 }
261
262 return fx*fy*fz;
263}
264
265
266
267
268Double_t AliTPCEfield::EvalField(Int_t ifun, Double_t x, Double_t y, Double_t z, Int_t type){
269 //
270 // Evaluate function ifun at position gx amd gy
271 // type == 0 - field
272 // == 1 - Ex
273 // == 2 - Ey
274 // == 3 - Ez
275 TMatrixD &mat = *fFitFunctions;
276 Int_t fid = TMath::Nint(mat(ifun,0));
277 Double_t ifx = (mat(ifun,1));
278 Double_t ify = (mat(ifun,2));
279 Double_t ifz = (mat(ifun,3));
280 //
281 if (type==0) return Field(fid,ifx,ify,ifz, x, y,z);
282 if (type>0) return FieldDn(fid,ifx,ify,ifz,type-1, x, y,z);
283 return 0;
284}
285
286Double_t AliTPCEfield::Eval(Double_t x, Double_t y, Double_t z, Int_t type){
287 //
288 // Evaluate function ifun at position gx amd gy
289 // type == 0 - field
290 // == 1 - Ex
291 // == 2 - Ey
292 // == 3 - Ez
293 Double_t value=0;
294 Double_t lx= 2.*(x-(fMin[0]+fMax[0])*0.5)/fScale;
295 Double_t ly= 2.*(y-(fMin[1]+fMax[1])*0.5)/fScale;
296 Double_t lz= 2.*(z-(fMin[2]+fMax[2])*0.5)/fScale;
297 //
298 Int_t nfun=fFitFunctions->GetNrows();
299 for (Int_t ifun=0; ifun<nfun; ifun++){
300 if (type==0) value+=(*fFitParam)[ifun]*EvalField(ifun,lx,ly,lz,type);
301 if (type>0) value+=2*(*fFitParam)[ifun]*EvalField(ifun,lx,ly,lz,type)/fScale;
302 }
303 return value;
304}
305
306Double_t AliTPCEfield::EvalS(Double_t x, Double_t y, Double_t z, Int_t type){
307 //
308 // static evaluation - possible to use it in the TF1
309 //
310 return fgInstance->Eval(x,y,z,type);
311}
312
313void AliTPCEfield::FitField(){
314 //
315 // Fit the e field
316 // Minimize chi2 residuals at the boundary points
317 // ?Tempoary sollution - integrals can be calculated analytically -
318 //
319 Int_t nfun=fFitFunctions->GetNrows();
320 Double_t *fun =new Double_t[nfun];
321 fFitter= new TLinearFitter(nfun, Form("hyp%d", nfun-1));
322 //
323 TTree * tree = GetTree("Boundary");
324 Int_t npoints = tree->GetEntries();
325 Int_t *indexes = new Int_t[npoints];
326 Double_t *rindex = new Double_t[npoints];
327 //
328
329 Double_t x=0, y=0, z=0, v=0;
330 tree->SetBranchAddress("x",&x);
331 tree->SetBranchAddress("y",&y);
332 tree->SetBranchAddress("z",&z);
333 tree->SetBranchAddress("v",&v);
334 TMatrixD valMatrix(npoints,4);
335 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
336 tree->GetEntry(ipoint);
337 valMatrix(ipoint,0)=2.*(x-(fMin[0]+fMax[0])*0.5)/fScale;
338 valMatrix(ipoint,1)=2.*(y-(fMin[1]+fMax[1])*0.5)/fScale;
339 valMatrix(ipoint,2)=2.*(z-(fMin[2]+fMax[2])*0.5)/fScale;
340 valMatrix(ipoint,3)=v;
341 rindex[ipoint]=gRandom->Rndm();
342 }
343 TMath::Sort(npoints,rindex,indexes);
344 for (Int_t jpoint=0; jpoint<npoints; jpoint++){
345 Int_t ipoint=indexes[jpoint];
346 //
347 Double_t lx= valMatrix(ipoint,0);
348 Double_t ly= valMatrix(ipoint,1);
349 Double_t lz= valMatrix(ipoint,2);
350 Double_t value = valMatrix(ipoint,3);
351 for (Int_t ifun=1; ifun<nfun; ifun++){
352 Double_t ffun=EvalField(ifun,lx,ly,lz,0);
353 fun[ifun-1]=ffun;
354 }
355 fFitter->AddPoint(fun,value,1);
356 }
357 fFitter->Eval();
358 fFitter->GetCovarianceMatrix(*fFitCovar);
359 fFitter->GetParameters(*fFitParam);
360
361 (*fWorkspace)<<"fitGlobal"<<
362 "covar.="<<fFitCovar<< // covariance matrix
363 "param.="<<fFitParam<< // fit parameters
364 "fun.="<<fFitFunctions<< // type of fit functions
365 "\n";
366
367 for (Int_t ifun=1; ifun<nfun; ifun++){
368
369 }
370
371
372}
373
374
375TMatrixD* AliTPCEfield::MakeCorrelation(TMatrixD &matrix){
376 //
377 //
378 //
379 Int_t nrows = matrix.GetNrows();
380 TMatrixD * mat = new TMatrixD(nrows,nrows);
381 for (Int_t irow=0; irow<nrows; irow++)
382 for (Int_t icol=0; icol<nrows; icol++){
383 (*mat)(irow,icol)= matrix(irow,icol)/TMath::Sqrt(matrix(irow,irow)*matrix(icol,icol));
384 }
385 return mat;
386}
387
388
389
390
391void AliTPCEfield::DumpField(Double_t gridSize, Double_t step){
392 //
393 //
394 //
395 Double_t stepSize=0.001*fScale/fMaxFreq;
396 //
397 for (Double_t x = fMin[0]+stepSize; x<=fMax[0]-stepSize; x+=gridSize){
398 for (Double_t y = fMin[1]+stepSize; y<=fMax[1]-stepSize; y+=gridSize)
399 for (Double_t z = fMin[2]; z<=fMax[2]; z+=gridSize){
400 //
401 //
402 Double_t v = Eval(x,y,z,0);
403 Double_t ex = Eval(x,y,z,1);
404 Double_t ey = Eval(x,y,z,2);
405 Double_t ez = Eval(x,y,z,3);
406 Double_t dexdx = (Eval(x,y,z,1)-Eval(x-stepSize,y,z,1))/stepSize; // numerical derivative
407 Double_t deydy = (Eval(x,y,z,2)-Eval(x,y-stepSize,z,2))/stepSize;
408 Double_t dezdz = (Eval(x,y,z,3)-Eval(x,y,z-stepSize,3))/stepSize;
409 (*fWorkspace)<<"dumpField"<<
410 "x="<<x<< // position
411 "y="<<y<<
412 "z="<<z<<
413 "v="<<v<< // potential
414 "ex="<<ex<< // Efield
415 "ey="<<ey<<
416 "ez="<<ez<<
417 "dexdx="<<dexdx<< //gradient of e field
418 "deydy="<<deydy<<
419 "dezdz="<<dezdz<<
420 "\n";
421 }
422 }
423 //
424 //
425 for (Double_t x = fMin[0]+stepSize; x<=fMax[0]-stepSize; x+=gridSize){
426 Double_t sumEx=0;
427 Double_t sumEy=0;
428 Double_t sumExEy=0;
429 for (Double_t y = fMin[1]+0.2; y<=fMax[1]-stepSize; y+=step)
430 for (Double_t z = fMin[2]; z<=fMax[2]; z+=gridSize){
431 //
432 //
433 Double_t v = Eval(x,y+step*0.5,z,0);
434 Double_t ex = Eval(x,y+step*0.5,z,1);
435 Double_t ey = Eval(x,y+step*0.5,z,2);
436 Double_t ez = Eval(x,y+step*0.5,z,3);
437 sumEx+=ex*step;
438 sumEy+=ey*step;
439 sumExEy+=(ex/ey)*step;
440 (*fWorkspace)<<"dumpDistortion"<<
441 "x="<<x<< // position
442 "y="<<y<<
443 "z="<<z<<
444 "v="<<v<< // potential
445 "ex="<<ex<< // Efield
446 "ey="<<ey<<
447 "ez="<<ez<<
448 //
449 "sEx="<<sumEx<< // field x integral
450 "sEy="<<sumEy<< // field y integral
451 "sExEy="<<sumExEy<< // tan integral
452 "\n";
453 }
454 }
455}
456
457
458void MakeTPC2DExample(AliTPCEfield *field){
459 //
460 /*
461 .L $ALICE_ROOT/TPC/AliTPCEfield.cxx++
462 AliTPCEfield *field = new AliTPCEfield("field",20, kTRUE,kTRUE);
463 MakeTPC2DExample(field)
464 field->FitField()
465 sqrt(field->fFitter.GetChisquare()/field->fFitter.GetNpoints())
466
467 TF2 f2("f2","AliTPCEfield::EvalS(x,y,0,0)",90,245,0,250);
468 f2->SetNpx(100); f2->SetNpy(100); f2->Draw("colz");
469
470 TF2 f2x("f2x","AliTPCEfield::EvalS(x,y,0,1)",90,240,0,240);
471 f2x->SetNpx(100); f2x->SetNpy(100); f2x->Draw("surf2");
472
473 TF2 f2y("f2y","AliTPCEfield::EvalS(x,y,0,2)",90,240,0,240);
474 f2y->SetNpx(100); f2y->SetNpy(100); f2y->Draw("surf2");
475
476 field->MakeCorrelation(*(field->fFitCovar)).Print()
477 Double_t index[100000];
478 for (Int_t i=0; i<field->fFitCovar->GetNrows(); i++) index[i]=i;
479 TGraph gr(field->fFitCovar->GetNrows(), index, field->fFitParam->GetMatrixArray());
480 gr->Draw("alp");
481
482 field->GetTree()->Draw("AliTPCEfield::EvalS(x,y,0,0):v");
483
484
485 TF2 f2xdy("f2xdy","AliTPCEfield::EvalS(x,y,0,1)/AliTPCEfield::EvalS(x,y,0,2)",90,240,0,240);
486 f2xdy->SetNpx(100); f2xdy->SetNpy(100); f2xdy->Draw("colz");
487
488 */
489
490 Double_t p0[4];
491 Double_t p1[4];
492 Double_t xmin=85, xmax=245;
493 Double_t ymin=0, ymax=250, deltaY=0.1*ymax;
494 Double_t vup=1;
495 Int_t npoints=1000;
496 field->SetRange(xmin, xmax,ymin,ymax,0,0);
497 // upper part
498 p0[0]=xmin; p0[1]=ymax+deltaY; p0[2]=0; p0[3]=vup;
499 p1[0]=xmax; p1[1]=ymax-deltaY; p1[2]=0; p1[3]=vup;
500 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],1,npoints);
501 //left
502 p0[0]=xmin; p0[1]=ymin+deltaY; p0[2]=0; p0[3]=0;
503 p1[0]=xmin; p1[1]=ymax+deltaY; p1[2]=0; p1[3]=vup;
504 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],2,npoints);
505 //right
506 p0[0]=xmax; p0[1]=ymin-deltaY; p0[2]=0; p0[3]=0;
507 p1[0]=xmax; p1[1]=ymax-deltaY; p1[2]=0; p1[3]=vup;
508 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],3,npoints);
509 //
510 //ROC
511 p0[0]=xmin; p0[1]=deltaY; p0[2]=0; p0[3]=-0;
512 p1[0]=xmax; p1[1]=-deltaY; p1[2]=0; p1[3]=0;
513 field->AddBoundaryLine(p0[0],p0[1], p0[2],p0[3],p1[0],p1[1], p1[2],p1[3],4,npoints);
514}
515
516
517