Script supersceeded by AliForwarddNdetaTask.C and
[u/mrichter/AliRoot.git] / TPC / AliTPCTempMap.cxx
CommitLineData
1209231c 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// //
19// TPC calibration class for temperature maps and tendencies //
20// (based on TPC Temperature Sensors and FiniteElement Simulation) //
21// //
22// Authors: Stefan Rossegger, Haavard Helstrup //
23// //
7ee86790 24// Note: Obvioulsy some changes by Marian, but when ??? //
25// //
1209231c 26///////////////////////////////////////////////////////////////////////////////
27
28#include "AliTPCSensorTempArray.h"
29#include "TLinearFitter.h"
30#include "TString.h"
31#include "TGraph2D.h"
f59cd9d0 32#include "TTimeStamp.h"
1209231c 33
34#include "AliTPCTempMap.h"
35
36
37ClassImp(AliTPCTempMap)
38
39 const char kStringFEsimulation[] = "FEsimulation.txt";
40
41//_____________________________________________________________________________
42AliTPCTempMap::AliTPCTempMap(AliTPCSensorTempArray *sensorDCS):
43 TNamed(),
f59cd9d0 44 fTempArray(0),
1209231c 45 fStringFEsimulation(kStringFEsimulation)
46{
47 //
48 // AliTPCTempMap default constructor
49 //
50
f59cd9d0 51 fTempArray = sensorDCS;
1209231c 52
53}
54
55//_____________________________________________________________________________
56AliTPCTempMap::AliTPCTempMap(const AliTPCTempMap &c):
57 TNamed(c),
f59cd9d0 58 fTempArray(c.fTempArray),
1209231c 59 fStringFEsimulation(c.fStringFEsimulation)
60{
61 //
62 // AliTPCTempMap copy constructor
63 //
64
65}
66
67//_____________________________________________________________________________
68AliTPCTempMap::~AliTPCTempMap()
69{
70 //
71 // AliTPCTempMap destructor
72 //
73
74}
75
76//_____________________________________________________________________________
77AliTPCTempMap &AliTPCTempMap::operator=(const AliTPCTempMap &c)
78{
79 //
80 // Assignment operator
81 //
82
83 if (this != &c) ((AliTPCTempMap &) c).Copy(*this);
84 return *this;
85
86}
87
88//_____________________________________________________________________________
89void AliTPCTempMap::Copy(TObject &c) const
90{
91 //
92 // Copy function
93 //
94
95 TObject::Copy(c);
96
97}
98
99//_____________________________________________________________________________
100
101Double_t AliTPCTempMap::GetTempGradientY(UInt_t timeSec, Int_t side){
102 //
103 // Extract Linear Vertical Temperature Gradient [K/cm] within the TPC on
104 // Shaft Side(A): 0
105 // Muon Side(C): 1
106 // Values based on TemperatureSensors within the TPC ( type: 3 (TPC) )
107 //
108 // FIXME: Also return residual-distribution, covariance Matrix
109 // or simply chi2 for validity check?
110 // -> better use GetLinearFitter - function in this case!
111
b479e253 112 TLinearFitter *fitter = new TLinearFitter(3,"x0++x1++x2");
1209231c 113 TVectorD param(3);
114 Int_t i = 0;
115
f59cd9d0 116 Int_t nsensors = fTempArray->NumSensors();
1209231c 117 for (Int_t isensor=0; isensor<nsensors; isensor++) { // loop over all sensors
f59cd9d0 118 AliTPCSensorTemp *entry = (AliTPCSensorTemp*)fTempArray->GetSensorNum(isensor);
1209231c 119
120 if (entry->GetType()==3 && entry->GetSide()==side) { // take SensorType:TPC
121 Double_t x[3];
122 x[0]=1;
123 x[1]=entry->GetX();
124 x[2]=entry->GetY();
f59cd9d0 125 Double_t y = fTempArray->GetValue(timeSec,isensor); // get temperature value
f4fe7830 126 if (IsOK(y)) fitter->AddPoint(x,y,1); // add values to LinearFitter
1209231c 127 i++;
128 }
129
130 }
b479e253 131 fitter->Eval();
132 fitter->GetParameters(param);
133
134 fitter->~TLinearFitter();
1209231c 135
136 return param[2]; // return vertical (Y) tempGradient in [K/cm]
137
138}
139
f1ea1647 140//_____________________________________________________________________________
f59cd9d0 141TLinearFitter *AliTPCTempMap::GetLinearFitter(Int_t type, Int_t side, TTimeStamp &stamp)
142{
143 //
144 // absolute time stamp used
145 // see AliTPCTempMap::GetLinearFitter(Int_t type, Int_t side, UInt_t timeSec) for details
146 //
148518a2 147 Int_t timeSec = stamp.GetSec()-fTempArray->GetStartTime().GetSec();
f59cd9d0 148 return GetLinearFitter(type,side,timeSec);
149}
150
1209231c 151//_____________________________________________________________________________
1209231c 152TLinearFitter *AliTPCTempMap::GetLinearFitter(Int_t type, Int_t side, UInt_t timeSec)
153{
154 //
155 // Creates a TlinearFitter object for the desired region of the TPC
156 // (via choosen type and side of TPC temperature sensors) at a given
157 // timeSec (in secounds) after start time
158 // type: 0 ... ReadOutChambers (ROC)
159 // 1 ... OuterContainmentVessel (OFC)
160 // 2 ... InnerContainmentVessel (IFC) + ThermalScreener (TS)
161 // 3 ... Within the TPC (DriftVolume) (TPC)
d0bd4fcc 162 // 4 ... Only InnerContainmentVessel (IFC)
1209231c 163 // side: Can be choosen for type 0 and 3 (otherwise it will be ignored in
164 // in order to get all temperature sensors of interest)
165 // 0 ... Shaft Side (A)
166 // 1 ... Muon Side (C)
167 //
168
169 TLinearFitter *fitter = new TLinearFitter(3);
170 Double_t *x = new Double_t[3];
171 Double_t y = 0;
f4fe7830 172 const Float_t kMaxDelta=0.5;
f59cd9d0 173
d0bd4fcc 174 if (type == 1 || type == 2 || type == 4) {
1209231c 175 fitter->SetFormula("x0++x1++TMath::Sin(x2)"); // returns Z,Y gradient
176 } else {
177 fitter->SetFormula("x0++x1++x2"); // returns X,Y gradient
178 }
179
180 Int_t i = 0;
f59cd9d0 181 Int_t nsensors = fTempArray->NumSensors();
182
183 Float_t temps[1000];
184 for (Int_t isensor=0; isensor<nsensors; isensor++) { // loop over all sensors
185 AliTPCSensorTemp *entry = (AliTPCSensorTemp*)fTempArray->GetSensorNum(isensor);
186 if (entry->GetType()==type && entry->GetSide()==side){
187 Float_t temperature= fTempArray->GetValue(timeSec,isensor); // get temperature value
f4fe7830 188 if (IsOK(temperature)) {temps[i]=temperature; i++;}
f59cd9d0 189 }
190 }
191 Float_t medianTemp = TMath::Median(i, temps);
192 if (i<3) return 0;
193 Float_t rmsTemp = TMath::RMS(i, temps);
194
195 i=0;
196
1209231c 197 for (Int_t isensor=0; isensor<nsensors; isensor++) { // loop over all sensors
f59cd9d0 198 AliTPCSensorTemp *entry = (AliTPCSensorTemp*)fTempArray->GetSensorNum(isensor);
1209231c 199
200 if (type==0 || type==3) { // 'side' information used
201 if (entry->GetType()==type && entry->GetSide()==side) {
202 x[0]=1;
203 x[1]=entry->GetX();
204 x[2]=entry->GetY();
f59cd9d0 205 y = fTempArray->GetValue(timeSec,isensor); // get temperature value
206 if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
f4fe7830 207 if (IsOK(y)) fitter->AddPoint(x,y,1); // add values to LinearFitter
1209231c 208 i++;
209 }
210 } else if (type==2) { // in case of IFC also usage of TS values
211 if ((entry->GetType()==2) || (entry->GetType()==5)) {
212 x[0]=1;
213 x[1]=entry->GetZ();
214 x[2]=entry->GetPhi();
f59cd9d0 215 y = fTempArray->GetValue(timeSec,isensor);
216 if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
f4fe7830 217 if (IsOK(y)) fitter->AddPoint(x,y,1);
1209231c 218 i++;
219 }
220 } else if (type==1){
221 if (entry->GetType()==type) {
222 x[0]=1;
223 x[1]=entry->GetZ();
224 x[2]=entry->GetPhi();
f59cd9d0 225 y = fTempArray->GetValue(timeSec,isensor);
226 if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
f4fe7830 227 if (IsOK(y)) fitter->AddPoint(x,y,1);
1209231c 228 i++;
229 }
d0bd4fcc 230 } else if (type==4) { // ONLY IFC
231 if (entry->GetType()==2) {
232 x[0]=1;
233 x[1]=entry->GetZ();
234 x[2]=entry->GetPhi();
f59cd9d0 235 y = fTempArray->GetValue(timeSec,isensor);
236 if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
f4fe7830 237 if (IsOK(y)) fitter->AddPoint(x,y,1);
d0bd4fcc 238 i++;
239 }
1209231c 240 }
241 }
f59cd9d0 242 fitter->Eval();
243 //fitter->EvalRobust(0.9); // Evaluates fitter
1209231c 244
245 delete [] x;
246
247 return fitter;
248
249 // returns TLinearFitter object where Chi2, Fitparameters and residuals can
250 // be extracted via usual memberfunctions
b479e253 251 // example: fitter->GetParameters(param)
1209231c 252 // In case of type IFC or OFC, the parameters are the gradients in
253 // Z and Y direction (see fitformula)
254 // Caution: Parameters are [K/cm] except Y at IFC,OFC ([K/radius])
255}
256
257//_____________________________________________________________________________
258
259TGraph2D *AliTPCTempMap::GetTempMapsViaSensors(Int_t type, Int_t side, UInt_t timeSec)
260{
261 //
262 // Creates a TGraph2D object for the desired region of the TPC
263 // (via choosen type and side of TPC temperature sensors) at a given
264 // timeSec (in secounds) after start time
265 // type: 0 ... ReadOutChambers (ROC)
266 // 1 ... OuterContainmentVessel (OFC)
267 // 2 ... InnerContainmentVessel (IFC) + ThermalScreener (TS)
268 // 3 ... Within the TPC (DriftVolume) (TPC)
269 // side: Can be choosen for type 0 and 3 (otherwise it will be ignored in
270 // in order to get all temperature sensors of interest)
55f06b51 271 // 0 ... A - side
272 // 1 ... C - side
1209231c 273 //
274
275 TGraph2D *graph2D = new TGraph2D();
276
277 Int_t i = 0;
278
f59cd9d0 279 Int_t nsensors = fTempArray->NumSensors();
1209231c 280 for (Int_t isensor=0; isensor<nsensors; isensor++) { // loop over all sensors
f59cd9d0 281 AliTPCSensorTemp *entry = (AliTPCSensorTemp*)fTempArray->GetSensorNum(isensor);
1209231c 282
283 Double_t x, y, z, r, phi, tempValue;
284 x = entry->GetX();
285 y = entry->GetY();
286 z = entry->GetZ();
287 r = entry->GetR();
288 phi = entry->GetPhi();
f59cd9d0 289 tempValue = fTempArray->GetValue(timeSec,isensor);
55f06b51 290 // printf("%d type %d: x=%lf y=%lf temp=%lf\n",isensor,entry->GetType(),x,y, tempValue);
1209231c 291 if (type==0 || type==3) { // 'side' information used
292 if (entry->GetType()==type && entry->GetSide()==side) {
293 graph2D->SetPoint(i,x,y,tempValue);
294 i++;
295 }
296 } else if (type==2) { // in case of IFC also usage of TS values
297 if (entry->GetType()==2 || entry->GetType()==5) {
298 graph2D->SetPoint(i,z,phi,tempValue);
299 i++;
300 }
301 } else if (type==1){
302 if (entry->GetType()==type) {
303 graph2D->SetPoint(i,z,phi,tempValue);
304 i++;
305 }
306 }
307 }
308
309 if (type==0 || type==3) {
310 graph2D->GetXaxis()->SetTitle("X[cm]");
311 graph2D->GetYaxis()->SetTitle("Y[cm]");
312 if (type==0 && side==0) {
55f06b51 313 graph2D->SetTitle("ROC A side");
1209231c 314 } else if (type==0 && side==1) {
55f06b51 315 graph2D->SetTitle("ROC C side");
1209231c 316 } else if (type==3 && side==0) {
55f06b51 317 graph2D->SetTitle("TPC A side (Inside the TPC)");
1209231c 318 } else if (type==3 && side==1) {
55f06b51 319 graph2D->SetTitle("TPC C side (Inside the TPC)");
1209231c 320 }
321 } else if (type==1 || type==2) {
322 graph2D->GetXaxis()->SetTitle("Z[cm]");
323 graph2D->GetYaxis()->SetTitle("Phi[RAD]");
324 if (type==1) {
325 graph2D->SetTitle("Outer Containment Vessel");
326 } else if (type==2) {
55f06b51 327 graph2D->SetTitle("Inner Containment Vessel");
1209231c 328 }
329 }
330
331 if (!graph2D->GetN()) {
332 printf("Returned TGraph2D is empty: check type and side values\n");
333 }
334
335 graph2D->GetXaxis()->SetLabelOffset(0.0);
336 graph2D->GetYaxis()->SetLabelOffset(0.005);
337 graph2D->GetZaxis()->SetLabelOffset(-0.04);
338
339
340 return graph2D; // returns TGgraph2D object
341
342}
343
344
345//_____________________________________________________________________________
346
347TGraph *AliTPCTempMap::MakeGraphGradient(Int_t axis, Int_t side, Int_t nPoints)
348{
349 //
350 // Make graph from start time to end time of TempGradient in axis direction
351 // axis: 0 ... horizontal Temperature Gradient (X)
352 // 1 ... vertical Temperature Gradient (Y)
353 // 2 ... longitudenal Temperature Gradient (Z) (side is ignored)
354 // z gradient value based on OFC temperature sensors
355 // Caution!: better z gradient values through difference between
356 // param[0] A- and param[0] C-side !
357 // side for X and Y gradient:
358 // 0 ... Shaft Side (A)
359 // 1 ... Muon Side (C)
360 //
361
362 TVectorD param(3);
363 TLinearFitter *fitter = new TLinearFitter(3);
364
f59cd9d0 365 UInt_t fStartTime = fTempArray->AliTPCSensorTempArray::GetStartTime();
366 UInt_t fEndTime = fTempArray->AliTPCSensorTempArray::GetEndTime();
1209231c 367
368 UInt_t stepTime = (fEndTime-fStartTime)/nPoints;
369
370 Double_t *x = new Double_t[nPoints];
371 Double_t *y = new Double_t[nPoints];
372 for (Int_t ip=0; ip<nPoints; ip++) {
373 x[ip] = fStartTime+ip*stepTime;
374 if (axis==2) {// Gradient in Z direction (based on OFC tempSensors)
375 fitter = GetLinearFitter(1, side, ip*stepTime);
376 } else {// Gradient in X or Y direction (based on TPC tempSensors)
377 fitter = GetLinearFitter(3, side, ip*stepTime);
378 }
379 fitter->GetParameters(param);
380 // multiplied by 500 since TempGradient is in [K/cm]
381 // (TPC diameter and length ~500cm)
382 if (axis==1) { // Y axis
383 y[ip] = param[2]*500;
384 } else { // X axis
385 y[ip] = param[1]*500;
386 }
387 }
388
389 TGraph *graph = new TGraph(nPoints,x,y);
390
391 fitter->~TLinearFitter();
392 delete [] x;
393 delete [] y;
394
395 graph->GetXaxis()->SetTimeDisplay(1);
396 graph->GetXaxis()->SetLabelOffset(0.02);
397 graph->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
398
399 return graph;
400}
401
f1ea1647 402
403//_____________________________________________________________________________
404Double_t AliTPCTempMap::GetTemperature(Double_t x, Double_t y, Double_t z, TTimeStamp &stamp)
405{
406 //
407 // absolute time stamp used
408 // see also Double_t AliTPCTempMap::GetTemperature(Double_t x, Double_t y, Double_t z, UInt_t timeSec) for details
409 //
410
411 Int_t timeSec = stamp.GetSec()-fTempArray->GetStartTime().GetSec();
412 return GetTemperature(x, y, z, timeSec);
413}
414
1209231c 415//_____________________________________________________________________________
416
417Double_t AliTPCTempMap::GetTemperature(Double_t x, Double_t y, Double_t z, UInt_t timeSec)
418{
419 //
f1ea1647 420 // Returns estimated Temperature at given position (x,y,z[cm]) at given time
1209231c 421 // (timeSec) after starttime
422 // Method: so far just a linear interpolation between Linar fits of
423 // the TPC temperature sensors
424 // FIXME: 'Educated Fit' through FiniteElement Simulation results!
425 // FIXXME: Return 0? if x,y,z out of range
426 //
427
428 TVectorD paramA(3), paramC(3);
429 TLinearFitter *fitterA = new TLinearFitter(3);
430 TLinearFitter *fitterC = new TLinearFitter(3);
431
432 fitterA = GetLinearFitter(3, 0, timeSec);
433 fitterA->GetParameters(paramA);
434 fitterC = GetLinearFitter(3, 1, timeSec);
435 fitterC->GetParameters(paramC);
436
437 Double_t fvalA = paramA[0]+paramA[1]*x+paramA[2]*y;
438 Double_t fvalC = paramC[0]+paramC[1]*x+paramC[2]*y;
439
440 Double_t k = (fvalA-fvalC)/(2*247);
441 Double_t tempValue = fvalC+(fvalA-fvalC)/2+k*z;
442
443 fitterA->~TLinearFitter();
444 fitterC->~TLinearFitter();
445
446 return tempValue;
f1ea1647 447
1209231c 448}
449
f4fe7830 450
451Bool_t AliTPCTempMap::IsOK(Float_t value){
452 //
7ee86790 453 // checks if value is within a certain range
f4fe7830 454 //
455 const Float_t kMinT=15;
456 const Float_t kMaxT=25;
457 return (value>kMinT && value<kMaxT);
458}