Test for Coverity
[u/mrichter/AliRoot.git] / TPC / AliTPCEfield.cxx
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
33 ClassImp(AliTPCEfield)
34
35
36 AliTPCEfield* AliTPCEfield::fgInstance=0;
37
38 AliTPCEfield::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
52 AliTPCEfield::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
69 void 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
125 AliTPCEfield::~AliTPCEfield() {
126   //
127   // Destructor
128   //
129   if (fWorkspace) delete fWorkspace;
130 }
131
132 void 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
145 void 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
173 TTree * AliTPCEfield::GetTree(const char * tname){
174   //  
175   //
176   //
177   return ((*fWorkspace)<<tname).GetTree();
178 }
179
180 Double_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
218 Double_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
268 Double_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
286 Double_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
306 Double_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
313 void 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
375 TMatrixD* 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
391 void 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
458 void 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