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