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