TPC/AliTPCExBEffective.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 // AliTPCExBEffective class                                                   //
  19 // Correct for the rest of ExB effect which are not covered by physical models
  20 //
  21 // Motivation:
  22 //   ExB correction:
  23 //      dr    =  c0* integral(Er/Ez) + c1* integral(Erphi/Ez)
  24 //      drphi = -c1* integral(Er/Ez) + c0* integral(Erphi/Ez)
  25 //   Where:
  26 //   wt = Bz*(k*vdrift/E)           ~ 0.3 at B=0.5 T
  27 //   c0 = 1/(1+T2*T2*wt*wt)
  28 //   c1 = T1*wt/(1+T1*T1*wt*wt)
  29 //
  30 // Residual integral(Er/Ez,Erphi/Ez) obtained comparing the B field 0 and B field +-0.5 T setting
  31 // minimizing track matching residuals
  32 // delta(Er/Ez) ~ sum[ poln(r) * polm(z) * cos(n,phi)]
  33 //
  34 ////////////////////////////////////////////////////////////////////////////
  35 #include "AliMagF.h"
  36 #include "TGeoGlobalMagField.h"
  37 #include "AliTPCcalibDB.h"
  38 #include "AliTPCParam.h"
  39 #include "AliLog.h"
  40
  41 #include "TMath.h"
  42 #include "AliTPCROC.h"
  43 #include "AliTPCExBEffective.h"
  44 ClassImp(AliTPCExBEffective)
  45
  46 AliTPCExBEffective::AliTPCExBEffective()
  47   : AliTPCCorrection("ExB_effective","ExB effective"),
  48     fC0(1.),fC1(0.),
  49     fPolynomA(0),
  50     fPolynomC(0),
  51     fPolynomValA(0),
  52     fPolynomValC(0)
  53 {
  54   //
  55   // default constructor
  56   //
  57 }
  58
  59 AliTPCExBEffective::~AliTPCExBEffective() {
  60   //
  61   // default destructor
  62   //
  63 }
  64
  65
  66
  67 void AliTPCExBEffective::Init() {
  68   //
  69   // Initialization funtion
  70   //
  71
  72   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  73   if (!magF) AliError("Magneticd field - not initialized");
  74   Double_t bzField = magF->SolenoidField()/10.; //field in T
  75   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
  76   if (!param) AliError("Parameters - not initialized");
  77   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
  78   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  79   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
  80   // Correction Terms for effective omegaTau; obtained by a laser calibration run
  81   SetOmegaTauT1T2(wt,fT1,fT2);
  82
  83
  84 }
  85
  86 void AliTPCExBEffective::Update(const TTimeStamp &/*timeStamp*/) {
  87   //
  88   // Update function
  89   //
  90   AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
  91   if (!magF) AliError("Magneticd field - not initialized");
  92   Double_t bzField = magF->SolenoidField()/10.; //field in T
  93   AliTPCParam *param= AliTPCcalibDB::Instance()->GetParameters();
  94   if (!param) AliError("Parameters - not initialized");
  95   Double_t vdrift = param->GetDriftV()/1000000.; // [cm/us]   // From dataBase: to be updated: per second (ideally)
  96   Double_t ezField = 400; // [V/cm]   // to be updated: never (hopefully)
  97   Double_t wt = -10.0 * (bzField*10) * vdrift / ezField ;
  98   // Correction Terms for effective omegaTau; obtained by a laser calibration run
  99   SetOmegaTauT1T2(wt,fT1,fT2);
 100
 101
 102 }
 103
 104
 105
 106 void AliTPCExBEffective::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 107   //
 108   // Calculates the correction due conical shape
 109   //
 110   if (!fPolynomA) return;
 111   AliTPCROC * calROC = AliTPCROC::Instance();
 112   const Double_t kRTPC0  =calROC->GetPadRowRadii(0,0);
 113   const Double_t kRTPC1  =calROC->GetPadRowRadii(36,calROC->GetNRows(36)-1);
 114   Float_t rmiddle=(kRTPC0+kRTPC1)/2.;
 115   //
 116   Double_t phi      = TMath::ATan2(x[1],x[0]);
 117   Double_t r        = TMath::Sqrt(x[1]*x[1]+x[0]*x[0]);
 118   Double_t driftN   = 1.-TMath::Abs(x[2])/calROC->GetZLength(0);  // drift from 0 to 1
 119   Double_t localxN  = 2*(r-rmiddle)/(kRTPC1-kRTPC0);         // normalize local x position
 120   //
 121   Double_t erez = 0;
 122   Double_t erphiez = 0;
 123   if (roc%36<18)  erez= GetSum(*fPolynomA, *fPolynomValA, localxN, driftN, phi,0);
 124   if (roc%36>=18) erez= GetSum(*fPolynomC, *fPolynomValC, localxN, driftN, phi,0);
 125   if (roc%36<18)  erphiez= GetSum(*fPolynomA, *fPolynomValA, localxN, driftN, phi,1);
 126   if (roc%36>=18) erphiez= GetSum(*fPolynomC, *fPolynomValC, localxN, driftN, phi,1);
 127
 128   Double_t dr    =   fC0 * erez + fC1 * erphiez;
 129   Double_t drphi =  -fC1 * erez + fC0 * erphiez;
 130
 131   // Calculate distorted position
 132   if ( r > 0.0 ) {
 133     r   =  r   + dr;
 134     phi =  phi + drphi/r;
 135   }
 136   // Calculate correction in cartesian coordinates
 137   dx[0] = r * TMath::Cos(phi) - x[0];
 138   dx[1] = r * TMath::Sin(phi) - x[1];
 139   dx[2] = 0.; // z distortion not implemented (1st order distortions)
 140
 141 }
 142
 143
 144
 145 Double_t AliTPCExBEffective::GetSum(const TMatrixD& mpol, const TMatrixD&mcoef, Double_t r, Double_t drift, Double_t phi, Int_t coord) const {
 146   //
 147   // Summation of the polynomials
 148   //
 149   Int_t npols=mpol.GetNrows();
 150   Double_t sum=0;
 151   for (Int_t ipol=0;ipol<npols; ipol++){
 152     Double_t pR = 1, pD=1, pPhi=1;
 153     Int_t icoord   = TMath::Nint(mpol(ipol,0));
 154     if (icoord!=coord) continue;
 155     Int_t npolR    = TMath::Nint(mpol(ipol,1));
 156     Int_t npolD    = TMath::Nint(mpol(ipol,2));
 157     Int_t npolPhi  = TMath::Nint(mpol(ipol,3));
 158     Double_t coef=mcoef(ipol,0);
 159     //
 160     for (Int_t ipolR=1; ipolR<=npolR; ipolR++) pR*=r;           // use simple polynoms
 161     for (Int_t ipolD=1; ipolD<=npolD; ipolD++) pD*=drift;       // use simple polynoms
 162     pPhi=TMath::Cos(npolPhi*phi);
 163     sum+= pR*pD*pPhi*coef;
 164   }
 165   return sum;
 166 }
 167
 168
 169 void AliTPCExBEffective::SetPolynoms(const TMatrixD *polA,const TMatrixD *polC){
 170   //
 171   // Set correction polynom - coefficients
 172   //
 173   fPolynomA = new TMatrixD(*polA);
 174   fPolynomC = new TMatrixD(*polC);
 175 }
 176
 177 void AliTPCExBEffective::SetCoeficients(const TMatrixD *valA,const TMatrixD *valC){
 178   //
 179   // Set correction polynom - coefficients
 180   //
 181   fPolynomValA = new TMatrixD(*valA);
 182   fPolynomValC = new TMatrixD(*valC);
 183 }
 184
 185
 186
 187
 188 void AliTPCExBEffective::Print(const Option_t* option) const {
 189   //
 190   // Print function to check the settings (e.g. the twist in the X direction)
 191   // option=="a" prints the C0 and C1 coefficents for calibration purposes
 192   //
 193
 194   TString opt = option; opt.ToLower();
 195   printf("%s\t%s\n",GetName(),GetTitle());
 196
 197   if (opt.Contains("a")) { // Print all details
 198     printf(" - T1: %1.4f, T2: %1.4f \n",fT1,fT2);
 199     printf(" - C0: %1.4f, C1: %1.4f \n",fC0,fC1);
 200     fPolynomValA->Print();
 201     fPolynomValC->Print();
 202   }
 203
 204
 205 }