[u/mrichter/AliRoot.git] / TPC / AliTPCTransformation.cxx

/*
  
 Class AliTPCtransformation:
 Should represent general non linear transformation. Currently tune for TPConly.
 To be used:
 1. Simulation-Digitization
 2. Reconstruction - AliTPCTransform
 3. Calibration/Alignment (KalmanFilter, Milipedde)
 4. Set of transformation to be stored/retrieved as OCDB entry 

 Base functionality:
 
 1. Double_t GetDeltaXYZ(Int_t coord, Int_t volID, Double_t param, Double_t x, Double_t y, Double_t z)
 Get correction - return the delta of coordinate coord dx or dy or dz for given volID for point at point (x,y,z)
 All coordinates are global

 2. The transformation should work only for given volIDs and detector IDs
    Currently Bitmask is used for filtering 


*/

#include <string.h>
#include "TRandom.h"
#include "TMath.h"
#include "TBits.h"
#include "TFormula.h"
#include "TF1.h"
#include "TLinearFitter.h"
#include "TFile.h"
#include "TObjString.h"

#include "TTreeStream.h"
#include "AliTrackPointArray.h"
#include "AliLog.h"
#include "AliTPCTransformation.h"

ClassImp(AliTPCTransformation)


AliTPCTransformation::GenFuncG    AliTPCTransformation::fgFormulas[10000];
TObjArray*  AliTPCTransformation::fgFormulasName = new TObjArray(10000);


void AliTPCTransformation::RegisterFormula(const char * name, GenFuncG formula){
  //
  // Add Formula to the list of formulas
  //
  Int_t last= fgFormulasName->GetEntries();
  fgFormulasName->AddLast(new TObjString(name));
  fgFormulas[last]=formula;
}

Int_t  AliTPCTransformation::BuildBasicFormulas(){
  //

  //
  //build list of basic TPC formulas - corrections
  //
  RegisterFormula("TPCscalingRPol",(GenFuncG)(AliTPCTransformation::TPCscalingRPol));
  RegisterFormula("TPCscalingRIFC",(GenFuncG)(AliTPCTransformation::TPCscalingRIFC));
  RegisterFormula("TPCscalingROFC",(GenFuncG)(AliTPCTransformation::TPCscalingROFC));
  //
  RegisterFormula("TPCscalingZDr",(GenFuncG)(AliTPCTransformation::TPCscalingZDr));
  RegisterFormula("TPCscalingPhiLocal",(GenFuncG)(AliTPCTransformation::TPCscalingPhiLocal));
  //
  // TPC Local X and Y misalignment 
  //
  RegisterFormula("TPClocaldLxdGX",(GenFuncG)(AliTPCTransformation::TPClocaldLxdGX));
  RegisterFormula("TPClocaldLxdGY",(GenFuncG)(AliTPCTransformation::TPClocaldLxdGY));
  RegisterFormula("TPClocaldLydGX",(GenFuncG)(AliTPCTransformation::TPClocaldLydGX));
  RegisterFormula("TPClocaldLydGY",(GenFuncG)(AliTPCTransformation::TPClocaldLydGY));
  //
  // Z offset
  //
  RegisterFormula("TPCDeltaZ",(GenFuncG)(AliTPCTransformation::TPCDeltaZ));
  RegisterFormula("TPCDeltaZMediumLong",(GenFuncG)(AliTPCTransformation::TPCDeltaZMediumLong));
  RegisterFormula("TPCTiltingZ",(GenFuncG)(AliTPCTransformation::TPCTiltingZ));
  return 0;
}

AliTPCTransformation::GenFuncG  AliTPCTransformation::FindFormula(const char * name){
  //
  // find formula - if registered
  //
  if (fgFormulasName->FindObject(name)==0) return 0;
  Int_t entries = fgFormulasName->GetEntries();
  for (Int_t i=0;i<entries;i++){
    if (strcmp(fgFormulasName->At(i)->GetName(), name)==0){
      return fgFormulas[i];
    }
  }
  return 0;
}

Double_t AliTPCTransformation::Eval(const char * name, const Double_t*x,const Double_t*par){
  //
  // Only for test purposes - very slow
  //
  GenFuncG fun = FindFormula(name);
  if (!fun) return 0;
  return fun(x,par);
}


AliTPCTransformation::AliTPCTransformation():
  TNamed(),
  fNameX(0),          // x formula name
  fNameY(0),          // y formula name
  fNameZ(0),          // z formula name 
  //
  fBitMask(0),        // bitmaps - transformation only for specified volID
  fCoordSystem(0),    // coord system of  output deltas 
  fParam(0),          // free parameter of transformation 
  fSigma(0),          // uncertainty of the parameter
  fSigma2Time(0),     // change of the error in time - (For kalman filter) 
  fFixedParam(0),     // fixed parameters of tranformation  
  fIsActive(kTRUE),   // swith - On/Off
  //
  fInit(kFALSE),      // initialization flag - set to kTRUE if corresponding formulas found
  fFormulaX(0),       // x formula - pointer to the function
  fFormulaY(0),       // y formula - pointer to the function
  fFormulaZ(0)       // z formula - pointer to the function
  //
{
  //
  // default constructor
  //
}


AliTPCTransformation::AliTPCTransformation(const char *name, TBits *mask, const char *fx, const char *fy, const char *fz, Int_t coordSystem):
  TNamed(name,name),
  fNameX(0),       // x formula name
  fNameY(0),       // y formula name
  fNameZ(0),       // z formula name
  fBitMask(mask),   // bitmaps - transformation only for specified volID
  fCoordSystem(coordSystem), // coordinate system of output deltas
  fParam(0),          // free parameter of transformation 
  fSigma(0),
  fSigma2Time(0),     // change of sigma in time
  fFixedParam(0),     // fixed parameters of tranformation  
  fIsActive(kTRUE),   // swith - On/Off
  //
  fInit(kFALSE),      // initialization flag - set to kTRUE if corresponding formulas found
  fFormulaX(0),       // x formula - pointer to the function
  fFormulaY(0),       // y formula - pointer to the function
  fFormulaZ(0)       // z formula - pointer to the function
{
  //
  // non default constructor
  //
  if (fx) fNameX= new TString(fx);
  if (fy) fNameY= new TString(fy);
  if (fz) fNameZ= new TString(fz);
  Init();
}

AliTPCTransformation::AliTPCTransformation(const AliTPCTransformation&trafo):
  TNamed(trafo),
  fNameX(0),          // x formula name
  fNameY(0),          // y formula name
  fNameZ(0),          // z formula name 
		     //
  fBitMask(0),        // bitmaps - transformation only for specified volID
  fCoordSystem(0),    // coord system of  output deltas 
  fParam(trafo.fParam),          // free parameter of transformation 
  fSigma(trafo.fSigma),          // uncertainty of the parameter
  fSigma2Time(trafo.fSigma2Time),     // change of the error in time - (For kalman filter) 
  fFixedParam(0),     // fixed parameters of tranformation  
  //
  fInit(kFALSE),      // initialization flag - set to kTRUE if corresponding formulas found
  fFormulaX(0),       // x formula - pointer to the function
  fFormulaY(0),       // y formula - pointer to the function
  fFormulaZ(0)       // z formula - pointer to the function
{
  if (trafo.fNameX) fNameX = new TString(*(trafo.fNameX)); 
  if (trafo.fNameY) fNameY = new TString(*(trafo.fNameY)); 
  if (trafo.fNameZ) fNameZ = new TString(*(trafo.fNameZ)); 
  if (trafo.fBitMask)  fBitMask = new TBits(*(trafo.fBitMask)); 
}

AliTPCTransformation::~AliTPCTransformation(){
  //
  // destructor
  //
  delete fNameX;
  delete fNameY;
  delete fNameZ;
  delete fBitMask;
  delete fFixedParam;
}
void AliTPCTransformation::SetParams(Double_t param, Double_t sigma, Double_t sigma2Time, TVectorD* fixedParams){
  //
  // Set parameters of transformation
  //
  fParam = param;
  fSigma = sigma;
  fSigma2Time = sigma2Time;
  if (fFixedParam) delete fFixedParam;
  fFixedParam = new TVectorD(*fixedParams);
  Init();
}


Bool_t AliTPCTransformation::Init(){
  //
  // associate formulas with pointer to the function
  //
  Bool_t isOK=kTRUE;
  if (fNameX) {
    fFormulaX=FindFormula(fNameX->Data());
    if (fFormulaX==0) isOK=kFALSE;
  }
  if (fNameY) {
    fFormulaY=FindFormula(fNameY->Data());
    if (fFormulaY==0) isOK=kFALSE;
  }
  if (fNameZ) {
    fFormulaZ=FindFormula(fNameZ->Data());
    if (!fFormulaZ) isOK=kFALSE;
  }
  return isOK;
}


TBits * AliTPCTransformation::BitsSide(Bool_t aside){
  //
  TBits * bits = new TBits(72);
  for (Int_t i=0; i<72;i++){
    if (i%36<18 && aside) (*bits)[i]=kTRUE;
    if (i%36<18 && (!aside)) (*bits)[i]=kFALSE;
    if (i%36>=18 && aside) (*bits)[i]=kFALSE;
    if (i%36>=18 && (!aside)) (*bits)[i]=kTRUE;
  }
  return bits;
}

TBits * AliTPCTransformation::BitsAll(){
  //
  //
  //
  TBits * bits = new TBits(72);
  for (Int_t i=0; i<72;i++){
    (*bits)[i]=kTRUE;
  }
  return bits;
}

Double_t AliTPCTransformation::GetDeltaXYZ(Int_t coord, Int_t volID, Double_t param, Double_t x, Double_t y, Double_t z){
  //
  //
  //
  if (!fIsActive) return 0;
  if (fBitMask && (!(*fBitMask)[volID])) return 0;
  Double_t xyz[5]={x,y,z, param,volID};
  if (fCoordSystem==0){
    // cartezian system
    if (coord==0 && fFormulaX) return fFormulaX(xyz,fFixedParam->GetMatrixArray()); 
    if (coord==1 && fFormulaY) return fFormulaY(xyz,fFixedParam->GetMatrixArray()); 
    if (coord==2 && fFormulaZ) return fFormulaZ(xyz,fFixedParam->GetMatrixArray());
  }
  if (fCoordSystem==1){  
    // cylindrical system
    if (coord==2) {
      if (fFormulaZ==0) return 0;
      return fFormulaZ(xyz,fFixedParam->GetMatrixArray());
    }
    Double_t rrphiz[3]={0,0,0};
    if (fFormulaX) rrphiz[0] = fFormulaX(xyz,fFixedParam->GetMatrixArray());
    if (fFormulaY) rrphiz[1] = fFormulaY(xyz,fFixedParam->GetMatrixArray());
    Double_t alpha = TMath::ATan2(y,x);
    Double_t ca    = TMath::Cos(alpha);
    Double_t sa    = TMath::Sin(alpha);
    if (coord==0) return ca*rrphiz[0]-sa*rrphiz[1];
    if (coord==1) return sa*rrphiz[0]+ca*rrphiz[1];
  }
  return 0;
}

Double_t  AliTPCTransformation::TPCscalingRPol(Double_t *xyz, Double_t * param){
  //
  // Scaling and shift of TPC radius
  // xyz[0..2] - global xyz of point 
  // xyz[3]    - scale parameter
  // param[0]  - radial scaling power
  // param[1]  - drift  scaling power
  // radius  from -1(at rInner)   to 1 (rOuter)
  // driftM  from -1(at 0 drift)  to 1 (250 cm drift)

  Double_t rInner=78.8;
  Double_t rOuter=258.0; 
  Double_t deltaR  = rOuter-rInner;  
  Double_t radius  = (TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])-rInner)*2./deltaR; 
  Double_t driftM  = (0.5 - TMath::Abs(xyz[2]/250.))*2.0;
  Double_t delta   = TMath::Power(radius,param[0])*TMath::Power(driftM,param[1]);
  return delta*xyz[3];
}


Double_t  AliTPCTransformation::TPCscalingZDr(Double_t *xyz, Double_t * param){
  //
  //
  // Scaling and shift of TPC radius
  // xyz[0..2] - global xyz of point 
  // xyz[3]    - scale parameter
  Double_t driftP  = TMath::Power(1. - TMath::Abs(xyz[2]/250.), param[0]);
  Double_t deltaZ  = (xyz[2]>0) ? -driftP : driftP;
  return deltaZ*xyz[3];
}

Double_t  AliTPCTransformation::TPCscalingPhiLocal(Double_t *xyz, Double_t * param){
  //
  //
  // Scaling if the local y -phi
  // xyz[0..2] - global xyz of point 
  // xyz[3]    - scale parameter
  // value = 1 for ful drift length and parameter 1
  Double_t alpha       = TMath::ATan2(xyz[1],xyz[0]);
  Double_t sector      = TMath::Nint(9*alpha/TMath::Pi()-0.5);
  Double_t localAlpha  = (alpha-(sector+0.5)*TMath::Pi()/9.);
  Double_t radius      = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])/250.;
  //
  Double_t deltaAlpha  = radius*TMath::Power(2.*9.*localAlpha/TMath::Pi(),param[0]);
  return deltaAlpha*xyz[3];
}


Double_t       AliTPCTransformation::TPCscalingRIFC(Double_t *xyz, Double_t * param){
  //
  // inner field cage r distorion - proportinal to 1 over distance to the IFC
  // param[0] - drift polynom order
  // distortion at first pad row - is normalized to 
  Double_t rInner=78.8;
  Double_t rFirst=85.2; 
  Double_t deltaR  = rFirst-rInner;
  Double_t ndistR  = (TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])-rInner)/deltaR;
  Double_t driftM  = (0.5 - TMath::Abs(xyz[2]/250.))*2.;
  Double_t value   = TMath::Power(driftM,param[0])/ndistR;
  return xyz[3]*value;
}

Double_t       AliTPCTransformation::TPCscalingROFC(Double_t *xyz, Double_t * param){
  //
  // outer field cage r distorion - proportinal to 1 over distance to the OFC
  // param[0] - drift polynom order
  // driftM   - from -1 to 1 
  //
  Double_t rLast=245.8;
  Double_t rOuter=258.0;  
  Double_t deltaR  = rOuter-rLast;
  Double_t ndistR  = (rOuter-TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]))/deltaR;
  Double_t driftM  = (0.5 - TMath::Abs(xyz[2]/250.))*2.;
  Double_t value   = TMath::Power(driftM,param[0])/ndistR;
  return xyz[3]*value;
}


//
// TPC sector local misalignment 
//
Double_t AliTPCTransformation:: TPClocaldLxdGX(Double_t *xyz, Double_t * /*param*/){
  //
  // xyz - [0..2] - position 
  //       [3]    - scale parameter
  //       [4]    - volID
  // return delta in global coordiante system
  //
  Int_t    sector = TMath::Nint(xyz[4]);
  Double_t alpha  = TMath::Pi()*(sector+0.5)/9;
  Double_t ca     = TMath::Cos(alpha);
  //  Double_t sa     = TMath::Sin(alpha);
  return ca*xyz[3];    
}

Double_t AliTPCTransformation::TPClocaldLxdGY(Double_t *xyz, Double_t * /*param*/){
  //
  // xyz - [0..2] - position 
  //       [3]    - scale parameter
  //       [4]    - volID
  // return delta in global coordiante system
  //
  Int_t    sector = TMath::Nint(xyz[4]);
  Double_t alpha  = TMath::Pi()*(sector+0.5)/9;
  //Double_t ca     = TMath::Cos(alpha);
  Double_t sa     = TMath::Sin(alpha);
  return   sa*xyz[3];  
}

Double_t AliTPCTransformation:: TPClocaldLydGX(Double_t *xyz, Double_t * /*param*/){
  //
  // xyz - [0..2] - position 
  //       [3]    - scale parameter
  //       [4]    - volID
  // return delta in global coordiante system
  //
  Int_t    sector = TMath::Nint(xyz[4]);
  Double_t alpha  = TMath::Pi()*(sector+0.5)/9;
  //Double_t ca     = TMath::Cos(alpha);
  Double_t sa     = TMath::Sin(alpha);
  return            -sa*xyz[3];  
}

Double_t AliTPCTransformation::TPClocaldLydGY(Double_t *xyz, Double_t * /*param*/){
  //
  // xyz - [0..2] - position 
  //       [3]    - scale parameter
  //       [4]    - volID
  // return delta in global coordiante system
  //
  Int_t    sector = TMath::Nint(xyz[4]);
  Double_t alpha  = TMath::Pi()*(sector+0.5)/9;
  Double_t ca     = TMath::Cos(alpha);
  //Double_t sa     = TMath::Sin(alpha);
  return   ca*xyz[3];  
}


Double_t        AliTPCTransformation::TPCDeltaZ(Double_t *xyz, Double_t * /*param*/){
  //
  // xyz - [0..2] - position 
  //        [3]    - scale parameter
  //        [4]    - volID
  // return delta in global coordiante system
  //
  Int_t    sector = TMath::Nint(xyz[4]);
  Double_t signZ  = (sector%36<18) ? 1: -1;  // drift direction
  return signZ*xyz[3];     // IROC shift
}

Double_t        AliTPCTransformation::TPCDeltaZMediumLong(Double_t *xyz, Double_t * /*param*/){
  //
  // xyz - [0..2] - position 
  //        [3]    - scale parameter
  //        [4]    - volID
  // return delta in global coordinate system 
  //
  Int_t    sector = TMath::Nint(xyz[4]);
  Double_t signZ  = (sector%36<18) ? 1: -1;  // drift direction
  if    (sector<36) return 0;     
  //
  Double_t radius  = (TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1])); 
  const Double_t radiusLong = 198.1;
  Double_t sign   = (radius<radiusLong) ? 1:-1;
  return xyz[3]*sign*signZ;
}

Double_t       AliTPCTransformation::TPCTiltingZ(Double_t *xyz, Double_t * param){
  // xyz - [0..2] - position 
  //        [3]    - scale parameter
  //        [4]    - volID
  // param[0]      - n for cos
  // param[1]      - n for sin
  // return delta in global coordinate system 
  Double_t radius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
  const Double_t rFirst=85.2; 
  const Double_t rLast =245.8;
  Double_t radiusC  = (rFirst+rLast)*0.5;
  Double_t deltaR = 2.0*(radius-radiusC)/(rLast-rFirst);
  Double_t alpha       = TMath::ATan2(xyz[1],xyz[0]);
  
  if (param[0]>0) deltaR *= TMath::Cos(param[0]*alpha);
  if (param[1]>0) deltaR *= TMath::Sin(param[1]*alpha);
  return deltaR*xyz[3];
}
Commit	Line	Data
94685752	1	/*
	2
	3	Class AliTPCtransformation:
	4	Should represent general non linear transformation. Currently tune for TPConly.
	5	To be used:
	6	1. Simulation-Digitization
	7	2. Reconstruction - AliTPCTransform
	8	3. Calibration/Alignment (KalmanFilter, Milipedde)
	9	4. Set of transformation to be stored/retrieved as OCDB entry
	10
	11	Base functionality:
	12
	13	1. Double_t GetDeltaXYZ(Int_t coord, Int_t volID, Double_t param, Double_t x, Double_t y, Double_t z)
	14	Get correction - return the delta of coordinate coord dx or dy or dz for given volID for point at point (x,y,z)
	15	All coordinates are global
	16
	17	2. The transformation should work only for given volIDs and detector IDs
	18	Currently Bitmask is used for filtering
	19
	20
	21
	22	*/
	23
	24	#include <string.h>
	25	#include "TRandom.h"
	26	#include "TMath.h"
	27	#include "TBits.h"
	28	#include "TFormula.h"
	29	#include "TF1.h"
	30	#include "TLinearFitter.h"
	31	#include "TFile.h"
	32	#include "TObjString.h"
	33
	34	#include "TTreeStream.h"
	35	#include "AliTrackPointArray.h"
	36	#include "AliLog.h"
	37	#include "AliTPCTransformation.h"
	38
	39	ClassImp(AliTPCTransformation)
	40
	41
	42	AliTPCTransformation::GenFuncG AliTPCTransformation::fgFormulas[10000];
	43	TObjArray* AliTPCTransformation::fgFormulasName = new TObjArray(10000);
	44
	45
	46	void AliTPCTransformation::RegisterFormula(const char * name, GenFuncG formula){
	47	//
	48	// Add Formula to the list of formulas
	49	//
	50	Int_t last= fgFormulasName->GetEntries();
	51	fgFormulasName->AddLast(new TObjString(name));
	52	fgFormulas[last]=formula;
	53	}
	54
	55	Int_t AliTPCTransformation::BuildBasicFormulas(){
	56	//
6d438146	57
	58	//
	59	//build list of basic TPC formulas - corrections
94685752	60	//
94685752	61	RegisterFormula("TPCscalingRPol",(GenFuncG)(AliTPCTransformation::TPCscalingRPol));
6d438146	62	RegisterFormula("TPCscalingRIFC",(GenFuncG)(AliTPCTransformation::TPCscalingRIFC));
	63	RegisterFormula("TPCscalingROFC",(GenFuncG)(AliTPCTransformation::TPCscalingROFC));
	64	//
94685752	65	RegisterFormula("TPCscalingZDr",(GenFuncG)(AliTPCTransformation::TPCscalingZDr));
94685752	66	RegisterFormula("TPCscalingPhiLocal",(GenFuncG)(AliTPCTransformation::TPCscalingPhiLocal));
0e9efcbe	67	//
	68	// TPC Local X and Y misalignment
	69	//
	70	RegisterFormula("TPClocaldLxdGX",(GenFuncG)(AliTPCTransformation::TPClocaldLxdGX));
	71	RegisterFormula("TPClocaldLxdGY",(GenFuncG)(AliTPCTransformation::TPClocaldLxdGY));
	72	RegisterFormula("TPClocaldLydGX",(GenFuncG)(AliTPCTransformation::TPClocaldLydGX));
	73	RegisterFormula("TPClocaldLydGY",(GenFuncG)(AliTPCTransformation::TPClocaldLydGY));
	74	//
	75	// Z offset
	76	//
	77	RegisterFormula("TPCDeltaZ",(GenFuncG)(AliTPCTransformation::TPCDeltaZ));
	78	RegisterFormula("TPCDeltaZMediumLong",(GenFuncG)(AliTPCTransformation::TPCDeltaZMediumLong));
	79	RegisterFormula("TPCTiltingZ",(GenFuncG)(AliTPCTransformation::TPCTiltingZ));
94685752	80	return 0;
	81	}
	82
	83	AliTPCTransformation::GenFuncG AliTPCTransformation::FindFormula(const char * name){
	84	//
	85	// find formula - if registered
	86	//
	87	if (fgFormulasName->FindObject(name)==0) return 0;
	88	Int_t entries = fgFormulasName->GetEntries();
	89	for (Int_t i=0;i<entries;i++){
	90	if (strcmp(fgFormulasName->At(i)->GetName(), name)==0){
	91	return fgFormulas[i];
	92	}
	93	}
	94	return 0;
	95	}
	96
	97	Double_t AliTPCTransformation::Eval(const char * name, const Double_tx,const Double_tpar){
	98	//
	99	// Only for test purposes - very slow
	100	//
	101	GenFuncG fun = FindFormula(name);
	102	if (!fun) return 0;
	103	return fun(x,par);
	104	}
	105
	106
	107	AliTPCTransformation::AliTPCTransformation():
	108	TNamed(),
	109	fNameX(0), // x formula name
	110	fNameY(0), // y formula name
	111	fNameZ(0), // z formula name
	112	//
	113	fBitMask(0), // bitmaps - transformation only for specified volID
	114	fCoordSystem(0), // coord system of output deltas
	115	fParam(0), // free parameter of transformation
	116	fSigma(0), // uncertainty of the parameter
6d438146	117	fSigma2Time(0), // change of the error in time - (For kalman filter)
94685752	118	fFixedParam(0), // fixed parameters of tranformation
0e9efcbe	119	fIsActive(kTRUE), // swith - On/Off
94685752	120	//
	121	fInit(kFALSE), // initialization flag - set to kTRUE if corresponding formulas found
	122	fFormulaX(0), // x formula - pointer to the function
	123	fFormulaY(0), // y formula - pointer to the function
	124	fFormulaZ(0) // z formula - pointer to the function
	125	//
	126	{
	127	//
	128	// default constructor
	129	//
	130	}
	131
	132
	133
6d438146	134	AliTPCTransformation::AliTPCTransformation(const char name, TBits mask, const char fx, const char fy, const char *fz, Int_t coordSystem):
94685752	135	TNamed(name,name),
	136	fNameX(0), // x formula name
	137	fNameY(0), // y formula name
	138	fNameZ(0), // z formula name
	139	fBitMask(mask), // bitmaps - transformation only for specified volID
	140	fCoordSystem(coordSystem), // coordinate system of output deltas
6d438146	141	fParam(0), // free parameter of transformation
	142	fSigma(0),
	143	fSigma2Time(0), // change of sigma in time
94685752	144	fFixedParam(0), // fixed parameters of tranformation
0e9efcbe	145	fIsActive(kTRUE), // swith - On/Off
94685752	146	//
	147	fInit(kFALSE), // initialization flag - set to kTRUE if corresponding formulas found
	148	fFormulaX(0), // x formula - pointer to the function
	149	fFormulaY(0), // y formula - pointer to the function
	150	fFormulaZ(0) // z formula - pointer to the function
	151	{
	152	//
	153	// non default constructor
	154	//
	155	if (fx) fNameX= new TString(fx);
	156	if (fy) fNameY= new TString(fy);
	157	if (fz) fNameZ= new TString(fz);
94685752	158	Init();
	159	}
	160
6d438146	161	AliTPCTransformation::AliTPCTransformation(const AliTPCTransformation&trafo):
	162	TNamed(trafo),
	163	fNameX(0), // x formula name
	164	fNameY(0), // y formula name
	165	fNameZ(0), // z formula name
	166	//
	167	fBitMask(0), // bitmaps - transformation only for specified volID
	168	fCoordSystem(0), // coord system of output deltas
	169	fParam(trafo.fParam), // free parameter of transformation
	170	fSigma(trafo.fSigma), // uncertainty of the parameter
	171	fSigma2Time(trafo.fSigma2Time), // change of the error in time - (For kalman filter)
	172	fFixedParam(0), // fixed parameters of tranformation
	173	//
	174	fInit(kFALSE), // initialization flag - set to kTRUE if corresponding formulas found
	175	fFormulaX(0), // x formula - pointer to the function
	176	fFormulaY(0), // y formula - pointer to the function
	177	fFormulaZ(0) // z formula - pointer to the function
	178	{
	179	if (trafo.fNameX) fNameX = new TString(*(trafo.fNameX));
	180	if (trafo.fNameY) fNameY = new TString(*(trafo.fNameY));
	181	if (trafo.fNameZ) fNameZ = new TString(*(trafo.fNameZ));
	182	if (trafo.fBitMask) fBitMask = new TBits(*(trafo.fBitMask));
	183	}
	184
	185	AliTPCTransformation::~AliTPCTransformation(){
	186	//
	187	// destructor
	188	//
	189	delete fNameX;
	190	delete fNameY;
	191	delete fNameZ;
	192	delete fBitMask;
	193	delete fFixedParam;
	194	}
	195	void AliTPCTransformation::SetParams(Double_t param, Double_t sigma, Double_t sigma2Time, TVectorD* fixedParams){
	196	//
	197	// Set parameters of transformation
	198	//
	199	fParam = param;
	200	fSigma = sigma;
	201	fSigma2Time = sigma2Time;
	202	if (fFixedParam) delete fFixedParam;
	203	fFixedParam = new TVectorD(*fixedParams);
	204	Init();
	205	}
94685752	206
	207
	208	Bool_t AliTPCTransformation::Init(){
	209	//
	210	// associate formulas with pointer to the function
	211	//
	212	Bool_t isOK=kTRUE;
	213	if (fNameX) {
	214	fFormulaX=FindFormula(fNameX->Data());
	215	if (fFormulaX==0) isOK=kFALSE;
	216	}
	217	if (fNameY) {
	218	fFormulaY=FindFormula(fNameY->Data());
	219	if (fFormulaY==0) isOK=kFALSE;
	220	}
	221	if (fNameZ) {
	222	fFormulaZ=FindFormula(fNameZ->Data());
	223	if (!fFormulaZ) isOK=kFALSE;
	224	}
	225	return isOK;
	226	}
	227
0e9efcbe	228
0e9efcbe	229
94685752	230	TBits * AliTPCTransformation::BitsSide(Bool_t aside){
	231	//
	232	TBits * bits = new TBits(72);
	233	for (Int_t i=0; i<72;i++){
	234	if (i%36<18 && aside) (*bits)[i]=kTRUE;
	235	if (i%36<18 && (!aside)) (*bits)[i]=kFALSE;
	236	if (i%36>=18 && aside) (*bits)[i]=kFALSE;
	237	if (i%36>=18 && (!aside)) (*bits)[i]=kTRUE;
	238	}
	239	return bits;
	240	}
	241
	242	TBits * AliTPCTransformation::BitsAll(){
	243	//
	244	//
	245	//
	246	TBits * bits = new TBits(72);
	247	for (Int_t i=0; i<72;i++){
	248	(*bits)[i]=kTRUE;
	249	}
	250	return bits;
	251	}
	252
	253	Double_t AliTPCTransformation::GetDeltaXYZ(Int_t coord, Int_t volID, Double_t param, Double_t x, Double_t y, Double_t z){
	254	//
	255	//
	256	//
0e9efcbe	257	if (!fIsActive) return 0;
94685752	258	if (fBitMask && (!(*fBitMask)[volID])) return 0;
0e9efcbe	259	Double_t xyz[5]={x,y,z, param,volID};
94685752	260	if (fCoordSystem==0){
	261	// cartezian system
	262	if (coord==0 && fFormulaX) return fFormulaX(xyz,fFixedParam->GetMatrixArray());
	263	if (coord==1 && fFormulaY) return fFormulaY(xyz,fFixedParam->GetMatrixArray());
0e9efcbe	264	if (coord==2 && fFormulaZ) return fFormulaZ(xyz,fFixedParam->GetMatrixArray());
94685752	265	}
	266	if (fCoordSystem==1){
	267	// cylindrical system
	268	if (coord==2) {
	269	if (fFormulaZ==0) return 0;
	270	return fFormulaZ(xyz,fFixedParam->GetMatrixArray());
	271	}
	272	Double_t rrphiz[3]={0,0,0};
	273	if (fFormulaX) rrphiz[0] = fFormulaX(xyz,fFixedParam->GetMatrixArray());
	274	if (fFormulaY) rrphiz[1] = fFormulaY(xyz,fFixedParam->GetMatrixArray());
	275	Double_t alpha = TMath::ATan2(y,x);
	276	Double_t ca = TMath::Cos(alpha);
	277	Double_t sa = TMath::Sin(alpha);
	278	if (coord==0) return carrphiz[0]-sarrphiz[1];
	279	if (coord==1) return sarrphiz[0]+carrphiz[1];
	280	}
	281	return 0;
	282	}
	283
	284	Double_t AliTPCTransformation::TPCscalingRPol(Double_t xyz, Double_t param){
	285	//
	286	// Scaling and shift of TPC radius
	287	// xyz[0..2] - global xyz of point
	288	// xyz[3] - scale parameter
6d438146	289	// param[0] - radial scaling power
	290	// param[1] - drift scaling power
	291	// radius from -1(at rInner) to 1 (rOuter)
	292	// driftM from -1(at 0 drift) to 1 (250 cm drift)
	293
	294	Double_t rInner=78.8;
	295	Double_t rOuter=258.0;
	296	Double_t deltaR = rOuter-rInner;
	297	Double_t radius = (TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1])-rInner)*2./deltaR;
	298	Double_t driftM = (0.5 - TMath::Abs(xyz[2]/250.))*2.0;
	299	Double_t delta = TMath::Power(radius,param[0])*TMath::Power(driftM,param[1]);
	300	return delta*xyz[3];
94685752	301	}
	302
	303
	304	Double_t AliTPCTransformation::TPCscalingZDr(Double_t xyz, Double_t param){
	305	//
	306	//
	307	// Scaling and shift of TPC radius
	308	// xyz[0..2] - global xyz of point
	309	// xyz[3] - scale parameter
	310	Double_t driftP = TMath::Power(1. - TMath::Abs(xyz[2]/250.), param[0]);
	311	Double_t deltaZ = (xyz[2]>0) ? -driftP : driftP;
	312	return deltaZ*xyz[3];
	313	}
	314
	315	Double_t AliTPCTransformation::TPCscalingPhiLocal(Double_t xyz, Double_t param){
	316	//
	317	//
	318	// Scaling if the local y -phi
	319	// xyz[0..2] - global xyz of point
	320	// xyz[3] - scale parameter
6d438146	321	// value = 1 for ful drift length and parameter 1
94685752	322	Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
	323	Double_t sector = TMath::Nint(9*alpha/TMath::Pi()-0.5);
	324	Double_t localAlpha = (alpha-(sector+0.5)*TMath::Pi()/9.);
	325	Double_t radius = TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1])/250.;
6d438146	326	//
6d438146	327	Double_t deltaAlpha = radiusTMath::Power(2.9.*localAlpha/TMath::Pi(),param[0]);
94685752	328	return deltaAlpha*xyz[3];
	329	}
	330
	331
6d438146	332	Double_t AliTPCTransformation::TPCscalingRIFC(Double_t xyz, Double_t param){
	333	//
	334	// inner field cage r distorion - proportinal to 1 over distance to the IFC
	335	// param[0] - drift polynom order
	336	// distortion at first pad row - is normalized to
	337	Double_t rInner=78.8;
	338	Double_t rFirst=85.2;
	339	Double_t deltaR = rFirst-rInner;
	340	Double_t ndistR = (TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1])-rInner)/deltaR;
	341	Double_t driftM = (0.5 - TMath::Abs(xyz[2]/250.))*2.;
	342	Double_t value = TMath::Power(driftM,param[0])/ndistR;
	343	return xyz[3]*value;
	344	}
	345
	346	Double_t AliTPCTransformation::TPCscalingROFC(Double_t xyz, Double_t param){
	347	//
	348	// outer field cage r distorion - proportinal to 1 over distance to the OFC
	349	// param[0] - drift polynom order
	350	// driftM - from -1 to 1
	351	//
	352	Double_t rLast=245.8;
	353	Double_t rOuter=258.0;
	354	Double_t deltaR = rOuter-rLast;
	355	Double_t ndistR = (rOuter-TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]))/deltaR;
	356	Double_t driftM = (0.5 - TMath::Abs(xyz[2]/250.))*2.;
	357	Double_t value = TMath::Power(driftM,param[0])/ndistR;
	358	return xyz[3]*value;
	359	}
0e9efcbe	360
	361
	362	//
	363	// TPC sector local misalignment
	364	//
	365	Double_t AliTPCTransformation:: TPClocaldLxdGX(Double_t xyz, Double_t /param/){
	366	//
	367	// xyz - [0..2] - position
	368	// [3] - scale parameter
	369	// [4] - volID
	370	// return delta in global coordiante system
	371	//
	372	Int_t sector = TMath::Nint(xyz[4]);
	373	Double_t alpha = TMath::Pi()*(sector+0.5)/9;
	374	Double_t ca = TMath::Cos(alpha);
	375	// Double_t sa = TMath::Sin(alpha);
	376	return ca*xyz[3];
	377	}
	378
	379	Double_t AliTPCTransformation::TPClocaldLxdGY(Double_t xyz, Double_t /param/){
	380	//
	381	// xyz - [0..2] - position
	382	// [3] - scale parameter
	383	// [4] - volID
	384	// return delta in global coordiante system
	385	//
	386	Int_t sector = TMath::Nint(xyz[4]);
	387	Double_t alpha = TMath::Pi()*(sector+0.5)/9;
	388	//Double_t ca = TMath::Cos(alpha);
	389	Double_t sa = TMath::Sin(alpha);
	390	return sa*xyz[3];
	391	}
	392
	393	Double_t AliTPCTransformation:: TPClocaldLydGX(Double_t xyz, Double_t /param/){
	394	//
	395	// xyz - [0..2] - position
	396	// [3] - scale parameter
	397	// [4] - volID
	398	// return delta in global coordiante system
	399	//
	400	Int_t sector = TMath::Nint(xyz[4]);
	401	Double_t alpha = TMath::Pi()*(sector+0.5)/9;
	402	//Double_t ca = TMath::Cos(alpha);
	403	Double_t sa = TMath::Sin(alpha);
	404	return -sa*xyz[3];
	405	}
	406
	407	Double_t AliTPCTransformation::TPClocaldLydGY(Double_t xyz, Double_t /param/){
	408	//
	409	// xyz - [0..2] - position
	410	// [3] - scale parameter
	411	// [4] - volID
	412	// return delta in global coordiante system
	413	//
	414	Int_t sector = TMath::Nint(xyz[4]);
	415	Double_t alpha = TMath::Pi()*(sector+0.5)/9;
	416	Double_t ca = TMath::Cos(alpha);
	417	//Double_t sa = TMath::Sin(alpha);
	418	return ca*xyz[3];
	419	}
	420
	421
	422	Double_t AliTPCTransformation::TPCDeltaZ(Double_t xyz, Double_t /param/){
	423	//
424	// xyz - [0..2] - position
425	// [3] - scale parameter
426	// [4] - volID
427	// return delta in global coordiante system
428	//
429	Int_t sector = TMath::Nint(xyz[4]);
430	Double_t signZ = (sector%36<18) ? 1: -1; // drift direction
431	return signZ*xyz[3]; // IROC shift
432	}
433
434	Double_t AliTPCTransformation::TPCDeltaZMediumLong(Double_t xyz, Double_t /param/){
435	//
436	// xyz - [0..2] - position
437	// [3] - scale parameter
438	// [4] - volID
439	// return delta in global coordinate system
440	//
441	Int_t sector = TMath::Nint(xyz[4]);
442	Double_t signZ = (sector%36<18) ? 1: -1; // drift direction
443	if (sector<36) return 0;
444	//
445	Double_t radius = (TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]));
446	const Double_t radiusLong = 198.1;
447	Double_t sign = (radius<radiusLong) ? 1:-1;
448	return xyz[3]signsignZ;
449	}
450
451	Double_t AliTPCTransformation::TPCTiltingZ(Double_t xyz, Double_t param){
452	// xyz - [0..2] - position
453	// [3] - scale parameter
454	// [4] - volID
455	// param[0] - n for cos
456	// param[1] - n for sin
457	// return delta in global coordinate system
458	Double_t radius = TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
459	const Double_t rFirst=85.2;
460	const Double_t rLast =245.8;
461	Double_t radiusC = (rFirst+rLast)*0.5;
462	Double_t deltaR = 2.0*(radius-radiusC)/(rLast-rFirst);
463	Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
464
465	if (param[0]>0) deltaR = TMath::Cos(param[0]alpha);
466	if (param[1]>0) deltaR = TMath::Sin(param[1]alpha);
467	return deltaR*xyz[3];
468	}
469