1ec91b0fe0b4d220d40dca4d500afcf1bcf2852b
[u/mrichter/AliRoot.git] / TPC / Base / AliTPCCalibVdrift.cxx
1 /**************************************************************************
2  * Copyright(c) 2006-07, 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 //                                                                           //
19 // Class describing the Vdrift dependencies on E,T,P and GasComposition      //
20 // Authors: Stefan Rossegger, Haavard Helstrup                               //
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include "TSystem.h"
25 #include "TObject.h"
26 #include "TMath.h"
27 #include "AliTPCTempMap.h"
28 #include "AliTPCSensorTempArray.h"
29
30 #include "AliTPCCalibVdrift.h"
31
32 ClassImp(AliTPCCalibVdrift)
33
34 namespace paramDefinitions {
35     
36   // Standard Conditions used as origin in the Magbolz simulations
37   // Dimesions E [kV/cm], T [K], P [TORR], Cco2 [%], Cn2 [%]
38   const Double_t kstdE = 400;
39   const Double_t kstdT = 293;
40   const Double_t kstdP = 744;
41   const Double_t kstdCco2 = 9.52;
42   const Double_t kstdCn2 = 4.76;
43   // Driftvelocity at Standardcontitions [cm/microSec]
44   const Double_t kstdVdrift = 2.57563;
45   
46   // Vdrift dependencies simulated with Magbolz [%(Vdrift)/[unit]]
47   const Double_t kdvdE = 0.24;
48   const Double_t kdvdT = 0.30;
49   const Double_t kdvdP = -0.13;
50   const Double_t kdvdCco2 = -6.60;
51   const Double_t kdvdCn2 = -1.74;
52   // 2nd order effect Taylor expansion
53   const Double_t kdvdE2nd = -0.00107628;
54   const Double_t kdvdT2nd = -0.00134441;
55   const Double_t kdvdP2nd = 0.000135325;
56   const Double_t kdvdCco22nd = 0.328761;
57   const Double_t kdvdCn22nd = 0.151605;
58
59   const Double_t torrTokPascal = 0.750061683;
60  
61   Double_t krho = 0.934246; // density of TPC-Gas [kg/m^3]
62                             // method of calculation: weighted average
63   Double_t kg = 9.81;
64
65   //
66   // Nominal value obtained from 2008 data
67   //
68   const Double_t kKelvin       =273.15; // degree to Kelvin
69   const Double_t kNominalTemp  =19.03;  // mean between A and C side  in degree
70   const Double_t kNominalPress =973.9;  // pressure sensor - in mbar- 
71                                         // calibDB->GetPressure(tstamp,irun,1)
72 }
73
74
75 using namespace paramDefinitions;
76
77 AliTPCCalibVdrift::AliTPCCalibVdrift():
78   TNamed(),
79   fSensTemp(0),
80   fSensPres(0),
81   fTempMap(0),
82   fSensGasComp(0),
83   fNominalTemp(0),    // nominal temperature in Kelvin
84   fNominalPress(0)    // nominal pressure    in mbar 
85 {
86   //
87   //  default constructor
88   //
89 }
90
91 AliTPCCalibVdrift::AliTPCCalibVdrift(AliTPCSensorTempArray *SensTemp, AliDCSSensor *SensPres, TObject *SensGasComp):
92   TNamed(),
93   fSensTemp(0),
94   fSensPres(0),
95   fTempMap(0),
96   fSensGasComp(0),
97   fNominalTemp(0),    // nominal temperature in Kelvin
98   fNominalPress(0)    // nominal pressure    in mbar 
99 {
100   //
101   //  Standard constructor
102   //
103
104   fSensTemp = SensTemp;
105   fSensPres = SensPres;
106   if (fSensTemp) {
107     fTempMap  = new AliTPCTempMap(fSensTemp);
108   } else {
109     fTempMap = 0;
110   }
111   fSensGasComp = SensGasComp;
112   fNominalTemp = kNominalTemp;
113   fNominalPress= kNominalPress;
114 }
115
116 //_____________________________________________________________________________
117 AliTPCCalibVdrift::AliTPCCalibVdrift(const AliTPCCalibVdrift& source) :
118   TNamed(source),
119   fSensTemp(source.fSensTemp),
120   fSensPres(source.fSensPres),
121   fTempMap(source.fTempMap),
122   fSensGasComp(source.fSensGasComp),
123   fNominalTemp(source.fNominalTemp),    // nominal temperature in Kelvin
124   fNominalPress(source.fNominalPress)    // nominal pressure    in mbar 
125
126 {
127   //
128   //  Copy constructor
129   //
130 }
131
132 //_____________________________________________________________________________
133 AliTPCCalibVdrift& AliTPCCalibVdrift::operator=(const AliTPCCalibVdrift& source){
134   //
135   // assignment operator
136   //
137   if (&source == this) return *this;
138   new (this) AliTPCCalibVdrift(source);
139   
140   return *this;  
141 }
142
143 //_____________________________________________________________________________
144 AliTPCCalibVdrift::~AliTPCCalibVdrift()
145 {
146   //
147   // AliTPCCalibVdrift destructor
148   // 
149
150 }
151
152 //_____________________________________________________________________________
153 Double_t AliTPCCalibVdrift::GetPTRelative(UInt_t absTimeSec, Int_t side){
154   //
155   // Get Relative difference of p/T for given time stamp
156   // absTimeSec - absolute time in secounds
157   // side: 0 - A side |  1 - C side
158   //
159
160   TTimeStamp tstamp(absTimeSec);
161
162   if (!fSensPres||!fSensTemp) return 0;
163   Double_t pressure = fSensPres->GetValue(tstamp);
164   TLinearFitter * fitter = fTempMap->GetLinearFitter(3,side,tstamp);
165   if (!fitter) return 0;
166   TVectorD vec;
167   fitter->GetParameters(vec);
168   delete fitter;
169   if (vec[0]<10) return 0;
170   //
171   //
172   //
173   Double_t  temperature = vec[0];  //vec[0] temeperature 
174   Double_t  tpnom       = (fNominalTemp+kKelvin)/(fNominalPress);
175   Double_t  tpmeasured  = (temperature+kKelvin)/(pressure);
176   Double_t  result      = (tpmeasured-tpnom)/tpnom;
177
178   return result;
179
180 }
181
182
183 //_____________________________________________________________________________
184 Double_t AliTPCCalibVdrift::VdriftLinearHyperplaneApprox(Double_t dE, Double_t dT, Double_t dP, Double_t dCco2, Double_t dCn2) 
185 {
186   //
187   // Returns approximated value for the driftvelocity change (in percent)
188   // based on a Hyperplane approximation (~ Taylorapproximation of 2nd order)
189   //
190
191   Double_t termE   = dE*kdvdE + TMath::Power(dE,2)*kdvdE2nd;
192   Double_t termT   = dT*kdvdT + TMath::Power(dT,2)*kdvdT2nd;
193   Double_t termP   = dP*kdvdP + TMath::Power(dP,2)*kdvdP2nd;
194   Double_t termCo2 = dCco2*kdvdCco2 + TMath::Power(dCco2,2)*kdvdCco22nd;
195   Double_t termN2  = dCn2*kdvdCn2 + TMath::Power(dCn2,2)*kdvdCn22nd;
196
197   Double_t vdChange = termE+termT+termP+termCo2+termN2;
198
199   return vdChange;
200
201 }
202
203 //_____________________________________________________________________________
204
205 Double_t AliTPCCalibVdrift::GetVdriftNominal() 
206 {
207   // returns nominal Driftvelocity at StandardConditions
208   return kstdVdrift;
209 }
210
211 //_____________________________________________________________________________
212
213 Double_t AliTPCCalibVdrift::GetVdriftChange(Double_t x, Double_t y, Double_t z, UInt_t absTimeSec)
214 {
215   // 
216   // Calculates Vdrift change in percent of Vdrift_nominal 
217   // (under nominal conditions) at x,y,z at absolute time (in sec)
218   //
219
220   TTimeStamp tstamp(absTimeSec);
221
222   // Get E-field Value --------------------------
223   Double_t dE = 0.23; // StandardOffset if CE is set to 100kV
224
225   // Get Temperature Value ----------------------  
226   AliTPCTempMap *tempMap = fTempMap;
227   Double_t dT = 0;
228   if (fTempMap) {
229     Double_t tempValue = tempMap->GetTemperature(x, y, z, tstamp);
230     dT = tempValue + 273.15 - kstdT;
231   }
232     
233   // Get Main Pressure Value ---------------------
234   Double_t dP = 0;
235   if (fSensPres==0) {
236     // Just the pressure drop over the TPC height
237     dP = - krho*kg*y/10000*torrTokPascal;
238   } else {
239     // pressure sensors plus additional 0.4mbar overpressure within the TPC
240     Double_t pressure = fSensPres->GetValue(tstamp) + 0.4; 
241     // calculate pressure drop according to height in TPC and transform to
242     // TORR (with simplified hydrostatic formula)
243     dP = (pressure - krho*kg*y/10000) * torrTokPascal - kstdP;
244   }
245
246   // Get GasComposition
247   // FIXME: include Goofy values for CO2 and N2 conzentration out of OCDB
248   //        Goofy not yet reliable ... 
249   Double_t dCco2 = 0;
250   Double_t dCn2 = 0;
251
252   // Calculate change in drift velocity in terms of Vdrift_nominal
253   Double_t vdChange = VdriftLinearHyperplaneApprox(dE, dT, dP, dCco2, dCn2); 
254   
255   return vdChange;
256     
257 }
258
259 //_____________________________________________________________________________
260
261 Double_t AliTPCCalibVdrift::GetMeanZVdriftChange(Double_t x, Double_t y, UInt_t absTimeSec)
262 {
263   // 
264   // Calculates Meanvalue in z direction of Vdrift change in percent 
265   // of Vdrift_nominal (under standard conditions) at position x,y,absTimeSec
266   // with help of 'nPopints' base points
267   //
268   
269   Int_t nPoints = 5;
270  
271   Double_t vdriftSum = 0;
272
273   for (Int_t i = 0; i<nPoints; i++) {
274     Double_t z = (Double_t)i/(nPoints-1)*500-250;
275     vdriftSum = vdriftSum + GetVdriftChange(x, y, z, absTimeSec);
276   }
277   
278   Double_t meanZVdrift = vdriftSum/nPoints;
279
280   return meanZVdrift;
281
282 }
283
284 //_____________________________________________________________________________
285
286 TGraph *AliTPCCalibVdrift::MakeGraphMeanZVdriftChange(Double_t x, Double_t y, Int_t nPoints)
287 {
288   //
289   // Make graph from start time to end time of Mean Drift Velocity in 
290   // Z direction at given x and y position
291   //
292
293   UInt_t startTime = fSensTemp->GetStartTime();
294   UInt_t endTime = fSensTemp->GetEndTime();
295   
296   UInt_t stepTime = (endTime - startTime)/nPoints;
297
298
299   Double_t *xvec = new Double_t[nPoints];
300   Double_t *yvec = new Double_t[nPoints];
301
302   for (Int_t ip=0; ip<nPoints; ip++) {
303     xvec[ip] = startTime+ip*stepTime;
304     yvec[ip] = GetMeanZVdriftChange(x, y, fSensTemp->GetStartTime().GetSec() + ip*stepTime);
305   }
306
307   TGraph *graph = new TGraph(nPoints,xvec,yvec);
308
309   delete [] xvec;
310   delete [] yvec;
311
312   graph->GetXaxis()->SetTimeDisplay(1);
313   graph->GetXaxis()->SetLabelOffset(0.02);
314   graph->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
315
316   return graph;
317 }