more secure string operations
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibVdrift.cxx
CommitLineData
1209231c 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
32ClassImp(AliTPCCalibVdrift)
33
da6c0bc9 34namespace paramDefinitions {
1209231c 35
f1ea1647 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;
1209231c 45
f1ea1647 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;
9430b11a 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)
1209231c 72}
73
f1ea1647 74
1209231c 75using namespace paramDefinitions;
76
da6c0bc9 77AliTPCCalibVdrift::AliTPCCalibVdrift(AliTPCSensorTempArray *SensTemp, AliDCSSensor *SensPres, TObject *SensGasComp):
1209231c 78 TNamed(),
79 fSensTemp(0),
80 fSensPres(0),
da6c0bc9 81 fTempMap(0),
9430b11a 82 fSensGasComp(0),
83 fNominalTemp(0), // nominal temperature in Kelvin
84 fNominalPress(0) // nominal pressure in mbar
1209231c 85{
86 //
87 // Standard constructor
88 //
89
90 fSensTemp = SensTemp;
91 fSensPres = SensPres;
f1ea1647 92 if (fSensTemp) {
93 fTempMap = new AliTPCTempMap(fSensTemp);
94 } else {
95 fTempMap = 0;
96 }
1209231c 97 fSensGasComp = SensGasComp;
9430b11a 98 fNominalTemp = kNominalTemp;
99 fNominalPress= kNominalPress;
1209231c 100}
101
f1ea1647 102//_____________________________________________________________________________
1209231c 103AliTPCCalibVdrift::AliTPCCalibVdrift(const AliTPCCalibVdrift& source) :
104 TNamed(source),
105 fSensTemp(source.fSensTemp),
106 fSensPres(source.fSensPres),
da6c0bc9 107 fTempMap(source.fTempMap),
9430b11a 108 fSensGasComp(source.fSensGasComp),
109 fNominalTemp(source.fNominalTemp), // nominal temperature in Kelvin
110 fNominalPress(source.fNominalPress) // nominal pressure in mbar
111
1209231c 112{
113 //
114 // Copy constructor
115 //
116}
117
118//_____________________________________________________________________________
1209231c 119AliTPCCalibVdrift& AliTPCCalibVdrift::operator=(const AliTPCCalibVdrift& source){
120 //
121 // assignment operator
122 //
123 if (&source == this) return *this;
124 new (this) AliTPCCalibVdrift(source);
125
126 return *this;
127}
128
129//_____________________________________________________________________________
130AliTPCCalibVdrift::~AliTPCCalibVdrift()
131{
132 //
133 // AliTPCCalibVdrift destructor
da6c0bc9 134 //
135
136}
137
f1ea1647 138//_____________________________________________________________________________
139Double_t AliTPCCalibVdrift::GetPTRelative(UInt_t absTimeSec, Int_t side){
1209231c 140 //
da6c0bc9 141 // Get Relative difference of p/T for given time stamp
f1ea1647 142 // absTimeSec - absolute time in secounds
143 // side: 0 - A side | 1 - C side
da6c0bc9 144 //
f1ea1647 145
146 TTimeStamp tstamp(absTimeSec);
147
148 if (!fSensPres||!fSensTemp) return 0;
149 Double_t pressure = fSensPres->GetValue(tstamp);
da6c0bc9 150 TLinearFitter * fitter = fTempMap->GetLinearFitter(3,side,tstamp);
151 if (!fitter) return 0;
152 TVectorD vec;
153 fitter->GetParameters(vec);
154 delete fitter;
155 if (vec[0]<10) return 0;
9430b11a 156 //
157 //
158 //
159 Double_t temperature = vec[0]; //vec[0] temeperature
160 Double_t tpnom = (fNominalTemp+kKelvin)/(fNominalPress);
161 Double_t tpmeasured = (temperature+kKelvin)/(pressure);
162 Double_t result = (tpmeasured-tpnom)/tpnom;
f1ea1647 163
9430b11a 164 return result;
f1ea1647 165
1209231c 166}
167
da6c0bc9 168
1209231c 169//_____________________________________________________________________________
170Double_t AliTPCCalibVdrift::VdriftLinearHyperplaneApprox(Double_t dE, Double_t dT, Double_t dP, Double_t dCco2, Double_t dCn2)
171{
172 //
f1ea1647 173 // Returns approximated value for the driftvelocity change (in percent)
174 // based on a Hyperplane approximation (~ Taylorapproximation of 2nd order)
1209231c 175 //
176
f1ea1647 177 Double_t termE = dE*kdvdE + TMath::Power(dE,2)*kdvdE2nd;
178 Double_t termT = dT*kdvdT + TMath::Power(dT,2)*kdvdT2nd;
179 Double_t termP = dP*kdvdP + TMath::Power(dP,2)*kdvdP2nd;
180 Double_t termCo2 = dCco2*kdvdCco2 + TMath::Power(dCco2,2)*kdvdCco22nd;
181 Double_t termN2 = dCn2*kdvdCn2 + TMath::Power(dCn2,2)*kdvdCn22nd;
182
183 Double_t vdChange = termE+termT+termP+termCo2+termN2;
184
185 return vdChange;
1209231c 186
187}
f1ea1647 188
1209231c 189//_____________________________________________________________________________
190
191Double_t AliTPCCalibVdrift::GetVdriftNominal()
192{
193 // returns nominal Driftvelocity at StandardConditions
194 return kstdVdrift;
195}
196
197//_____________________________________________________________________________
198
f1ea1647 199Double_t AliTPCCalibVdrift::GetVdriftChange(Double_t x, Double_t y, Double_t z, UInt_t absTimeSec)
1209231c 200{
201 //
202 // Calculates Vdrift change in percent of Vdrift_nominal
f1ea1647 203 // (under nominal conditions) at x,y,z at absolute time (in sec)
1209231c 204 //
205
f1ea1647 206 TTimeStamp tstamp(absTimeSec);
207
1209231c 208 // Get E-field Value --------------------------
f1ea1647 209 Double_t dE = 0.23; // StandardOffset if CE is set to 100kV
1209231c 210
211 // Get Temperature Value ----------------------
da6c0bc9 212 AliTPCTempMap *tempMap = fTempMap;
f1ea1647 213 Double_t dT = 0;
214 if (fTempMap) {
215 Double_t tempValue = tempMap->GetTemperature(x, y, z, tstamp);
216 dT = tempValue + 273.15 - kstdT;
217 }
218
1209231c 219 // Get Main Pressure Value ---------------------
f1ea1647 220 Double_t dP = 0;
221 if (fSensPres==0) {
222 // Just the pressure drop over the TPC height
223 dP = - krho*kg*y/10000*torrTokPascal;
224 } else {
225 // pressure sensors plus additional 0.4mbar overpressure within the TPC
226 Double_t pressure = fSensPres->GetValue(tstamp) + 0.4;
227 // calculate pressure drop according to height in TPC and transform to
228 // TORR (with simplified hydrostatic formula)
229 dP = (pressure - krho*kg*y/10000) * torrTokPascal - kstdP;
230 }
231
1209231c 232 // Get GasComposition
f1ea1647 233 // FIXME: include Goofy values for CO2 and N2 conzentration out of OCDB
234 // Goofy not yet reliable ...
1209231c 235 Double_t dCco2 = 0;
236 Double_t dCn2 = 0;
237
238 // Calculate change in drift velocity in terms of Vdrift_nominal
f1ea1647 239 Double_t vdChange = VdriftLinearHyperplaneApprox(dE, dT, dP, dCco2, dCn2);
1209231c 240
f1ea1647 241 return vdChange;
242
1209231c 243}
244
245//_____________________________________________________________________________
246
f1ea1647 247Double_t AliTPCCalibVdrift::GetMeanZVdriftChange(Double_t x, Double_t y, UInt_t absTimeSec)
1209231c 248{
249 //
250 // Calculates Meanvalue in z direction of Vdrift change in percent
f1ea1647 251 // of Vdrift_nominal (under standard conditions) at position x,y,absTimeSec
1209231c 252 // with help of 'nPopints' base points
253 //
254
255 Int_t nPoints = 5;
256
2aaa30ef 257 Double_t vdriftSum = 0;
1209231c 258
259 for (Int_t i = 0; i<nPoints; i++) {
260 Double_t z = (Double_t)i/(nPoints-1)*500-250;
f1ea1647 261 vdriftSum = vdriftSum + GetVdriftChange(x, y, z, absTimeSec);
1209231c 262 }
263
2aaa30ef 264 Double_t meanZVdrift = vdriftSum/nPoints;
1209231c 265
2aaa30ef 266 return meanZVdrift;
1209231c 267
268}
269
270//_____________________________________________________________________________
271
272TGraph *AliTPCCalibVdrift::MakeGraphMeanZVdriftChange(Double_t x, Double_t y, Int_t nPoints)
273{
274 //
275 // Make graph from start time to end time of Mean Drift Velocity in
276 // Z direction at given x and y position
277 //
278
2aaa30ef 279 UInt_t startTime = fSensTemp->GetStartTime();
280 UInt_t endTime = fSensTemp->GetEndTime();
1209231c 281
2aaa30ef 282 UInt_t stepTime = (endTime - startTime)/nPoints;
1209231c 283
284
285 Double_t *xvec = new Double_t[nPoints];
286 Double_t *yvec = new Double_t[nPoints];
287
288 for (Int_t ip=0; ip<nPoints; ip++) {
2aaa30ef 289 xvec[ip] = startTime+ip*stepTime;
f1ea1647 290 yvec[ip] = GetMeanZVdriftChange(x, y, fSensTemp->GetStartTime().GetSec() + ip*stepTime);
1209231c 291 }
292
293 TGraph *graph = new TGraph(nPoints,xvec,yvec);
294
295 delete [] xvec;
296 delete [] yvec;
297
298 graph->GetXaxis()->SetTimeDisplay(1);
299 graph->GetXaxis()->SetLabelOffset(0.02);
300 graph->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
301
302 return graph;
303}