coverity fix
[u/mrichter/AliRoot.git] / TPC / AliTPCExBEffective.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 //                                                                        //
18 // AliTPCExBEffective class                                                   //
19 // Correct for the rest of ExB effect which are not covered by physical models
20 //
21 // Motivation:
22 //   ExB correction: 
23 //      dr    =  c0* integral(Er/Ez) + c1* integral(Erphi/Ez)
24 //      drphi = -c1* integral(Er/Ez) + c0* integral(Erphi/Ez)
25 //   Where:      
26 //   wt = Bz*(k*vdrift/E)           ~ 0.3 at B=0.5 T 
27 //   c0 = 1/(1+T2*T2*wt*wt) 
28 //   c1 = T1*wt/(1+T1*T1*wt*wt)
29 //   
30 // Residual integral(Er/Ez,Erphi/Ez) obtained comparing the B field 0 and B field +-0.5 T setting
31 // minimizing track matching residuals 
32 // delta(Er/Ez) ~ sum[ poln(r) * polm(z) * cos(n,phi)] 
33 //  
34 ////////////////////////////////////////////////////////////////////////////
35 #include "AliMagF.h"
36 #include "TGeoGlobalMagField.h"
37 #include "AliTPCcalibDB.h"
38 #include "AliTPCParam.h"
39 #include "AliLog.h"
40
41 #include "TMath.h"
42 #include "AliTPCROC.h"
43 #include "AliTPCExBEffective.h"
44 ClassImp(AliTPCExBEffective)
45
46 AliTPCExBEffective::AliTPCExBEffective()
47   : AliTPCCorrection("ExB_effective","ExB effective"),
48     fC0(1.),fC1(0.), 
49     fPolynomA(0),
50     fPolynomC(0),
51     fPolynomValA(0),
52     fPolynomValC(0)
53 {
54   //
55   // default constructor
56   //
57 }
58
59 AliTPCExBEffective::~AliTPCExBEffective() {
60   //
61   // default destructor
62   //
63 }
64
65
66
67 void AliTPCExBEffective::Init() {
68   //
69   // Initialization funtion
70   //
71   
72   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
73   if (!magF) AliError("Magneticd field - not initialized");
74   Double_t bzField = magF->SolenoidField()/10.; //field in T
75   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
76   if (!param) AliError("Parameters - not initialized");
77   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
78   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
79   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
80   // Correction Terms for effective omegaTau; obtained by a laser calibration run
81   SetOmegaTauT1T2(wt,fT1,fT2);
82
83
84 }
85
86 void AliTPCExBEffective::Update(const TTimeStamp &/*timeStamp*/) {
87   //
88   // Update function 
89   //
90   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
91   if (!magF) AliError("Magneticd field - not initialized");
92   Double_t bzField = magF->SolenoidField()/10.; //field in T
93   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
94   if (!param) AliError("Parameters - not initialized");
95   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
96   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
97   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ; 
98   // Correction Terms for effective omegaTau; obtained by a laser calibration run
99   SetOmegaTauT1T2(wt,fT1,fT2);
100
101
102 }
103
104
105
106 void AliTPCExBEffective::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
107   //
108   // Calculates the correction due conical shape
109   //   
110   if (!fPolynomA) return;
111   AliTPCROC * calROC = AliTPCROC::Instance();
112   const Double_t kRTPC0  =calROC->GetPadRowRadii(0,0);
113   const Double_t kRTPC1  =calROC->GetPadRowRadii(36,calROC->GetNRows(36)-1);
114   Float_t rmiddle=(kRTPC0+kRTPC1)/2.;
115   //
116   Double_t phi      = TMath::ATan2(x[1],x[0]);
117   Double_t r        = TMath::Sqrt(x[1]*x[1]+x[0]*x[0]);
118   Double_t driftN   = 1.-TMath::Abs(x[2])/calROC->GetZLength(0);  // drift from 0 to 1
119   Double_t localxN  = 2*(r-rmiddle)/(kRTPC1-kRTPC0);         // normalize local x position
120   //
121   Double_t erez = 0;
122   Double_t erphiez = 0;
123   if (roc%36<18)  erez= GetSum(*fPolynomA, *fPolynomValA, localxN, driftN, phi,0);
124   if (roc%36>=18) erez= GetSum(*fPolynomC, *fPolynomValC, localxN, driftN, phi,0);
125   if (roc%36<18)  erphiez= GetSum(*fPolynomA, *fPolynomValA, localxN, driftN, phi,1);
126   if (roc%36>=18) erphiez= GetSum(*fPolynomC, *fPolynomValC, localxN, driftN, phi,1);
127
128   Double_t dr    =   fC0 * erez + fC1 * erphiez;
129   Double_t drphi =  -fC1 * erez + fC0 * erphiez;
130   
131   // Calculate distorted position
132   if ( r > 0.0 ) {
133     r   =  r   + dr;
134     phi =  phi + drphi/r;
135   }
136   // Calculate correction in cartesian coordinates
137   dx[0] = r * TMath::Cos(phi) - x[0];
138   dx[1] = r * TMath::Sin(phi) - x[1];
139   dx[2] = 0.; // z distortion not implemented (1st order distortions)
140
141 }
142
143
144
145 Double_t AliTPCExBEffective::GetSum(const TMatrixD& mpol, const TMatrixD&mcoef, Double_t r, Double_t drift, Double_t phi, Int_t coord) const {
146   //
147   // Summation of the polynomials
148   //
149   Int_t npols=mpol.GetNrows();
150   Double_t sum=0;
151   for (Int_t ipol=0;ipol<npols; ipol++){
152     Double_t pR = 1, pD=1, pPhi=1;
153     Int_t icoord   = TMath::Nint(mpol(ipol,0));
154     if (icoord!=coord) continue;
155     Int_t npolR    = TMath::Nint(mpol(ipol,1));
156     Int_t npolD    = TMath::Nint(mpol(ipol,2));
157     Int_t npolPhi  = TMath::Nint(mpol(ipol,3));
158     Double_t coef=mcoef(ipol,0);
159     //
160     for (Int_t ipolR=1; ipolR<=npolR; ipolR++) pR*=r;           // use simple polynoms
161     for (Int_t ipolD=1; ipolD<=npolD; ipolD++) pD*=drift;       // use simple polynoms
162     pPhi=TMath::Cos(npolPhi*phi);
163     sum+= pR*pD*pPhi*coef;
164   }
165   return sum;
166 }
167  
168
169 void AliTPCExBEffective::SetPolynoms(const TMatrixD *polA,const TMatrixD *polC){
170   //
171   // Set correction polynom - coefficients
172   //
173   fPolynomA = new TMatrixD(*polA);
174   fPolynomC = new TMatrixD(*polC);
175 }
176
177 void AliTPCExBEffective::SetCoeficients(const TMatrixD *valA,const TMatrixD *valC){
178   //
179   // Set correction polynom - coefficients
180   //
181   fPolynomValA = new TMatrixD(*valA);
182   fPolynomValC = new TMatrixD(*valC);
183 }
184
185
186
187
188 void AliTPCExBEffective::Print(const Option_t* option) const {
189   //
190   // Print function to check the settings (e.g. the twist in the X direction)
191   // option=="a" prints the C0 and C1 coefficents for calibration purposes
192   //
193
194   TString opt = option; opt.ToLower();
195   printf("%s\t%s\n",GetName(),GetTitle());
196   
197   if (opt.Contains("a")) { // Print all details
198     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
199     printf(" - C0: %1.4f, C1: %1.4f \n",fC0,fC1);
200     fPolynomValA->Print();
201     fPolynomValC->Print();
202   }    
203  
204  
205 }