more secure string operations
[u/mrichter/AliRoot.git] / TPC / AliTPCTempMap.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 //                                                                           //
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 //                                                                           //
24 //  Note: Obvioulsy some changes by Marian, but when ???                     //
25 //                                                                           //
26 ///////////////////////////////////////////////////////////////////////////////
27
28 #include "AliTPCSensorTempArray.h"
29 #include "TLinearFitter.h"
30 #include "TString.h"
31 #include "TGraph2D.h"
32 #include "TTimeStamp.h"
33
34 #include "AliTPCTempMap.h"
35
36
37 ClassImp(AliTPCTempMap)
38   
39   const char kStringFEsimulation[] = "FEsimulation.txt";
40
41 //_____________________________________________________________________________
42 AliTPCTempMap::AliTPCTempMap(AliTPCSensorTempArray *sensorDCS):
43   TNamed(),
44   fTempArray(0),
45   fStringFEsimulation(kStringFEsimulation)
46 {
47   //
48   // AliTPCTempMap default constructor
49   //
50
51   fTempArray = sensorDCS;
52
53 }
54
55 //_____________________________________________________________________________
56 AliTPCTempMap::AliTPCTempMap(const AliTPCTempMap &c):
57   TNamed(c),
58   fTempArray(c.fTempArray),
59   fStringFEsimulation(c.fStringFEsimulation)
60 {
61   //
62   // AliTPCTempMap copy constructor
63   //
64
65 }
66
67 //_____________________________________________________________________________
68 AliTPCTempMap::~AliTPCTempMap()
69 {
70   //
71   // AliTPCTempMap destructor
72   //
73   
74 }
75
76 //_____________________________________________________________________________
77 AliTPCTempMap &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 //_____________________________________________________________________________
89 void AliTPCTempMap::Copy(TObject &c) const
90 {
91   //
92   // Copy function
93   //
94   
95   TObject::Copy(c);
96   
97 }
98
99 //_____________________________________________________________________________
100
101 Double_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   
112  TLinearFitter *fitter = new TLinearFitter(3,"x0++x1++x2");
113  TVectorD param(3);
114  Int_t i = 0;
115
116  Int_t nsensors = fTempArray->NumSensors();
117  for (Int_t isensor=0; isensor<nsensors; isensor++) { // loop over all sensors
118    AliTPCSensorTemp *entry = (AliTPCSensorTemp*)fTempArray->GetSensorNum(isensor);
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();    
125      Double_t y = fTempArray->GetValue(timeSec,isensor); // get temperature value
126      if (IsOK(y)) fitter->AddPoint(x,y,1); // add values to LinearFitter
127      i++;
128    }
129
130  }  
131  fitter->Eval();
132  fitter->GetParameters(param);
133
134  fitter->~TLinearFitter();
135
136  return param[2]; // return vertical (Y) tempGradient in [K/cm]
137   
138 }
139
140 //_____________________________________________________________________________
141 TLinearFitter *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   //
147   Int_t timeSec = stamp.GetSec()-fTempArray->GetStartTime().GetSec();
148   return GetLinearFitter(type,side,timeSec);
149 }
150
151 //_____________________________________________________________________________
152 TLinearFitter *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)
162   //       4 ... Only InnerContainmentVessel (IFC) 
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;
172   const Float_t kMaxDelta=0.5;
173   
174   if (type == 1 || type == 2 || type == 4) {
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;
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
188       if (IsOK(temperature)) {temps[i]=temperature; i++;}
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   
197   for (Int_t isensor=0; isensor<nsensors; isensor++) { // loop over all sensors
198     AliTPCSensorTemp *entry = (AliTPCSensorTemp*)fTempArray->GetSensorNum(isensor);
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();    
205         y = fTempArray->GetValue(timeSec,isensor); // get temperature value
206         if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
207         if (IsOK(y)) fitter->AddPoint(x,y,1); // add values to LinearFitter
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();    
215         y = fTempArray->GetValue(timeSec,isensor);
216         if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
217         if (IsOK(y)) fitter->AddPoint(x,y,1); 
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();    
225         y = fTempArray->GetValue(timeSec,isensor);
226         if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
227         if (IsOK(y)) fitter->AddPoint(x,y,1);
228         i++;    
229       }
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();    
235         y = fTempArray->GetValue(timeSec,isensor);
236         if (TMath::Abs(y-medianTemp)>kMaxDelta+4.*rmsTemp) continue;
237         if (IsOK(y)) fitter->AddPoint(x,y,1); 
238         i++;
239       }
240     }
241   }  
242   fitter->Eval();
243   //fitter->EvalRobust(0.9); // Evaluates fitter
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
251   // example: fitter->GetParameters(param)
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
259 TGraph2D *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)
271   //       0 ... A - side
272   //       1 ... C - side
273   // 
274
275   TGraph2D *graph2D = new TGraph2D();
276
277   Int_t i = 0;
278   
279   Int_t nsensors = fTempArray->NumSensors();
280   for (Int_t isensor=0; isensor<nsensors; isensor++) { // loop over all sensors
281     AliTPCSensorTemp *entry = (AliTPCSensorTemp*)fTempArray->GetSensorNum(isensor);
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();
289     tempValue = fTempArray->GetValue(timeSec,isensor);
290     //    printf("%d type %d: x=%lf y=%lf temp=%lf\n",isensor,entry->GetType(),x,y, tempValue);
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) {
313       graph2D->SetTitle("ROC A side");
314     } else if (type==0 && side==1) {
315       graph2D->SetTitle("ROC C side");
316     } else if (type==3 && side==0) {
317       graph2D->SetTitle("TPC A side (Inside the TPC)");
318     } else if (type==3 && side==1) {
319       graph2D->SetTitle("TPC C side (Inside the TPC)");
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) {
327       graph2D->SetTitle("Inner Containment Vessel");
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
347 TGraph *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
365   UInt_t fStartTime = fTempArray->AliTPCSensorTempArray::GetStartTime();
366   UInt_t fEndTime = fTempArray->AliTPCSensorTempArray::GetEndTime();
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
402
403 //_____________________________________________________________________________
404 Double_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
415 //_____________________________________________________________________________
416
417 Double_t AliTPCTempMap::GetTemperature(Double_t x, Double_t y, Double_t z, UInt_t timeSec)
418 {  
419   //
420   // Returns estimated Temperature at given position (x,y,z[cm]) at given time 
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;
447
448 }
449
450
451 Bool_t  AliTPCTempMap::IsOK(Float_t value){
452   //
453   // checks if value is within a certain range
454   //
455   const Float_t kMinT=15;
456   const Float_t kMaxT=25;
457   return (value>kMinT && value<kMaxT);
458 }