Adding simple example to load default debug streamer
[u/mrichter/AliRoot.git] / TPC / AliTPCExBFirst.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 // This implementation AliTPCExB is using an aproximate calculation of the ExB
18 // effect. Therefore the drift ODE is Taylor expanded and only the first
19 // order part is taken.
20 //
21 // 
22 // The ExB correction map is stored in the calib DB
23 // To test current version:
24 /*
25   char *storage = "local://OCDBres"
26   Int_t RunNumber=0;
27   AliCDBManager::Instance()->SetDefaultStorage(storage);
28   AliCDBManager::Instance()->SetRun(RunNumber) 
29   AliTPCExBFirst * exb = AliTPCcalibDB::Instance()->GetExB();
30
31   //  exb->TestExB("exb.root");
32   // TFile f("exb.root");
33   //positions->Draw("drphi");
34 */
35
36
37
38 #include "TMath.h"
39 #include "AliFieldMap.h"
40 #include "AliMagF.h"
41 #include "TTreeStream.h"
42 #include "AliTPCExBFirst.h"
43
44 ClassImp(AliTPCExBFirst)
45
46 const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
47 const Double_t AliTPCExBFirst::fgkDriftField=40.e3;
48
49 AliTPCExBFirst::AliTPCExBFirst()
50   : fDriftVelocity(0),
51     fkNX(0),fkNY(0),fkNZ(0),
52     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
53     fkZMin(-250.),fkZMax(250.),
54     fkNMean(0),
55     fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
56   //
57   // purely for I/O
58   //
59   SetInstance(this);
60 }
61
62 AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField,
63                                Double_t driftVelocity,
64                                Int_t nx,Int_t ny,Int_t nz)
65   : fDriftVelocity(driftVelocity),
66     fkNX(nx),fkNY(ny),fkNZ(nz),
67     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
68     fkZMin(-250.),fkZMax(250.),
69     fkNMean(0),
70     fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
71   //
72   // The constructor. One has to supply a magnetic field and an (initial)
73   // drift velocity. Since some kind of lookuptable is created the
74   // number of its meshpoints can be supplied.
75   //
76   ConstructCommon(0,bField);
77   SetInstance(this);
78 }
79
80 AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
81                                Double_t driftVelocity) 
82   : fDriftVelocity(driftVelocity),
83     fkNX(0),fkNY(0),fkNZ(0),
84     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
85     fkZMin(-250.),fkZMax(250.),
86     fkNMean(0),
87     fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
88   //
89   // The constructor. One has to supply a field map and an (initial)
90   // drift velocity.
91   //
92   SetInstance(this);
93   fkXMin=bFieldMap->Xmin()
94     -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
95     *bFieldMap->DelX();
96   fkXMax=bFieldMap->Xmax()
97     -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
98     *bFieldMap->DelX();
99   fkYMin=bFieldMap->Ymin()
100     -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
101     *bFieldMap->DelY();
102   fkYMax=bFieldMap->Ymax()
103     -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
104     *bFieldMap->DelY();
105   fkZMin=bFieldMap->Zmin()
106     -TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
107     *bFieldMap->DelZ();
108   fkZMax=bFieldMap->Zmax()
109     -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
110     *bFieldMap->DelZ();
111
112   fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
113   fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
114   fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
115
116   ConstructCommon(bFieldMap,0);
117 }
118
119 AliTPCExBFirst::~AliTPCExBFirst() { 
120   //
121   // destruct the poor object.
122   //
123   delete[] fkMeanBx;
124   delete[] fkMeanBy;
125 }
126
127 void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
128   //
129   // correct for the distortion
130   //
131   Double_t bx,by;
132   GetMeanFields(position[0],position[1],position[2],&bx,&by);
133   if (position[2]>0.) {
134     Double_t bxe,bye;
135     GetMeanFields(position[0],position[1],250.,&bxe,&bye);
136     if (position[2]!=250.) {
137       bx=(500.*bxe-(position[2]+250.)*bx)/(250.-position[2]);
138       by=(500.*bye-(position[2]+250.)*by)/(250.-position[2]);
139     }
140     else {
141       bx=bxe;
142       by=bye;
143     }
144   }
145   
146   Double_t mu=fDriftVelocity/fgkDriftField;
147   Double_t wt=mu*fkMeanBz;
148   
149   corrected[0]=mu*(wt*bx-by)/(1.+wt*wt);
150   corrected[1]=mu*(wt*by+bx)/(1.+wt*wt);
151   
152   if (position[2]>0.) {
153     corrected[0]*=(250.-position[2]);
154     corrected[1]*=(250.-position[2]);
155   }
156   else {
157     corrected[0]*=(-250.-position[2]);
158     corrected[1]*=(-250.-position[2]);
159   }
160
161   corrected[0]=position[0]-corrected[0];
162   corrected[1]=position[1]-corrected[1];
163   corrected[2]=position[2];
164 }
165
166 void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
167   //
168   // well, as the name sais...
169   //
170   TTreeSRedirector ts(fileName);
171   Double_t x[3];
172   for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
173     for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
174       for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
175         Double_t d[3];
176         Correct(x,d);
177         Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
178         Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
179         Double_t dr=r-rd;
180         Double_t phi=TMath::ATan2(x[0],x[1]);
181         Double_t phid=TMath::ATan2(d[0],d[1]);
182         Double_t dphi=phi-phid;
183         if (dphi<0.) dphi+=TMath::TwoPi();
184         if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
185         Double_t drphi=r*dphi;
186         Double_t dx=x[0]-d[0];
187         Double_t dy=x[1]-d[1];
188         Double_t dz=x[2]-d[2];
189         ts<<"positions"
190           <<"x0="<<x[0]
191           <<"x1="<<x[1]
192           <<"x2="<<x[2]
193           <<"dx="<<dx
194           <<"dy="<<dy
195           <<"dz="<<dz
196           <<"r="<<r
197           <<"phi="<<phi
198           <<"dr="<<dr
199           <<"drphi="<<drphi
200           <<"\n";
201       }
202 }
203
204 void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap,
205                                      const AliMagF *bField) {
206   //
207   // THIS IS PRIVATE! (a helper for the constructor)
208   //
209   fkNMean=fkNX*fkNY*fkNZ;
210   fkMeanBx=new Double_t[fkNMean];
211   fkMeanBy=new Double_t[fkNMean];
212
213   Double_t x[3];
214   Double_t nBz=0;
215   fkMeanBz=0.;
216   for (int i=0;i<fkNX;++i) {
217     x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1);
218     for (int j=0;j<fkNY;++j) {
219       x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1);
220       Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
221       Double_t bx=0.,by=0.;
222       for (int k=0;k<fkNZ;++k) {
223         x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
224         Float_t b[3];
225         // the x is not const in the Field function...
226         Float_t xt[3];
227         for (int l=0;l<3;++l) xt[l]=x[l];
228         // that happens due to the lack of a sophisticated class design:
229         if (bFieldMap!=0)
230           bFieldMap->Field(xt,b);
231         else 
232           bField->Field(xt,b);
233         bx+=b[0]/10.;
234         by+=b[1]/10.;
235         fkMeanBx[(k*fkNY+j)*fkNX+i]=bx;
236         fkMeanBy[(k*fkNY+j)*fkNX+i]=by;
237         if (90.<=r&&r<=250.) {
238           fkMeanBz+=b[2]/10.;
239           ++nBz;
240         }
241       }
242     }
243   }
244   fkMeanBz/=nBz;
245 }
246
247 void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
248                                    Double_t *Bx,Double_t *By) const {
249   //
250   // THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
251   //
252   Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin);
253   Int_t xi1=static_cast<Int_t>(x);
254   xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
255   Int_t xi2=xi1+1;
256   Double_t dx=(x-xi1);
257   Double_t dx1=(xi2-x);
258
259   Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin);
260   Int_t yi1=static_cast<Int_t>(y);
261   yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
262   Int_t yi2=yi1+1;
263   Double_t dy=(y-yi1);
264   Double_t dy1=(yi2-y);
265   
266   Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin);
267   Int_t zi1=static_cast<Int_t>(z);
268   zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
269   Int_t zi2=zi1+1;
270   Double_t dz=(z-zi1);
271
272   double s0x=fkMeanBx[yi1*fkNX+xi1]*dx1*dy1
273             +fkMeanBx[yi2*fkNX+xi1]*dx1*dy
274             +fkMeanBx[yi1*fkNX+xi2]*dx *dy
275             +fkMeanBx[yi2*fkNX+xi2]*dx *dy1;
276   double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1
277             +fkMeanBy[yi2*fkNX+xi1]*dx1*dy
278             +fkMeanBy[yi1*fkNX+xi2]*dx *dy
279             +fkMeanBy[yi2*fkNX+xi2]*dx *dy1;
280   Int_t zi0=zi1-1;
281   double snmx,snmy;
282   if (zi0>=0) {
283     snmx=fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
284         +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
285         +fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
286         +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
287     snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
288         +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
289         +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
290         +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
291   }
292   else
293     snmx=snmy=0.;
294   double snx=fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
295             +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
296             +fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
297             +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
298   double sny=fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
299             +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
300             +fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
301             +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
302   double snpx=fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
303              +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
304              +fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
305              +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
306   double snpy=fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
307              +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
308              +fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
309              +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
310   *Bx=0.5*(((snpx-2.*snx+snmx)*dz+2.*(snx-snmx))*dz+snx-s0x+snmx);
311   *By=0.5*(((snpy-2.*sny+snmy)*dz+2.*(sny-snmy))*dz+sny-s0y+snmy);
312   //TODO: make this nice
313   if (TMath::Abs(z)>0.001) {
314     *Bx/=z;
315     *By/=z;
316   }
317 }