HLT/TPCLib/transform/AliHLTTPCSpline2D3D.cxx

   1 //**************************************************************************
   2 //* This file is property of and copyright by the ALICE HLT Project        *
   3 //* ALICE Experiment at CERN, All rights reserved.                         *
   4 //*                                                                        *
   5 //* Primary Authors: Sergey Gorbunov <sergey.gorbunov@cern.ch>             *
   6 //*                  for The ALICE HLT Project.                            *
   7 //*                                                                        *
   8 //* Permission to use, copy, modify and distribute this software and its   *
   9 //* documentation strictly for non-commercial purposes is hereby granted   *
  10 //* without fee, provided that the above copyright notice appears in all   *
  11 //* copies and that both the copyright notice and this permission notice   *
  12 //* appear in the supporting documentation. The authors make no claims     *
  13 //* about the suitability of this software for any purpose. It is          *
  14 //* provided "as is" without express or implied warranty.                  *
  15 //**************************************************************************
  16
  17 /** @file   AliHLTTPCSpline2D3D.cxx
  18     @author Sergey Gorbubnov
  19     @date
  20     @brief
  21 */
  22
  23
  24 #include "AliHLTTPCFastTransform.h"
  25 #include "AliTPCTransform.h"
  26 #include "AliTPCParam.h"
  27 #include "AliTPCcalibDB.h"
  28 #include "Vc/Vc"
  29
  30 #include <iostream>
  31 #include <iomanip>
  32
  33 using namespace std;
  34
  35
  36
  37 void AliHLTTPCSpline2D3D::Init(Float_t minA,Float_t  maxA, Int_t  nBinsA, Float_t  minB,Float_t  maxB, Int_t  nBinsB)
  38 {
  39   //
  40   // Initialisation
  41   //
  42
  43   if( maxA<= minA ) maxA = minA+1;
  44   if( maxB<= minB ) maxB = minB+1;
  45   if( nBinsA <4 ) nBinsA = 4;
  46   if( nBinsB <4 ) nBinsB = 4;
  47
  48   fNA = nBinsA;
  49   fNB = nBinsB;
  50   fN = fNA*fNB;
  51   fMinA = minA;
  52   fMinB = minB;
  53
  54   fStepA = (maxA-minA)/(nBinsA-1);
  55   fStepB = (maxB-minB)/(nBinsB-1);
  56   fScaleA = 1./fStepA;
  57   fScaleB = 1./fStepB;
  58
  59   Vc::free( fXYZ );
  60   fXYZ = Vc::malloc< float, Vc::AlignOnCacheline>( 4*fN );
  61   memset ( fXYZ, 0, fN*4*sizeof(Float_t) );
  62 }
  63
  64 AliHLTTPCSpline2D3D::~AliHLTTPCSpline2D3D()
  65 {
  66   Vc::free( fXYZ );
  67 }
  68
  69 void AliHLTTPCSpline2D3D::Consolidate()
  70 {
  71   //
  72   // Consolidate the map
  73   //
  74
  75   Float_t *m = fXYZ;
  76   for( int iXYZ=0; iXYZ<3; iXYZ++ ){
  77     for( int iA=0; iA<fNA; iA++){
  78       {
  79         int i0 = 4*iA*fNB + iXYZ;
  80         int i1 = i0+4;
  81         int i2 = i0+4*2;
  82         int i3 = i0+4*3;
  83         m[i0] = 0.5*( m[i3]+m[i0]+3*(m[i1]-m[i2]) );
  84       }
  85       {
  86         int i0 = 4*(iA*fNB+fNB-4) + iXYZ;
  87         int i1 = i0+4;
  88         int i2 = i0+4*2;
  89         int i3 = i0+4*3;
  90         m[i3] = 0.5*( m[i0]+m[i3]+3*(m[i2]-m[i1]) );
  91       }
  92     }
  93     for( int iB=0; iB<fNB; iB++){
  94       {
  95         int i0 = 4*(0*fNB +iB) +iXYZ;
  96         int i1 = i0+4*fNB;
  97         int i2 = i1+4*fNB;
  98         int i3 = i2+4*fNB;
  99         m[i0] = 0.5*( m[i3]+m[i0]+3*(m[i1]-m[i2]) );
 100       }
 101       {
 102         int i0 = 4*((fNA-4)*fNB +iB) +iXYZ;
 103         int i1 = i0+4*fNB;
 104         int i2 = i1+4*fNB;
 105         int i3 = i2+4*fNB;
 106         m[i3] = 0.5*( m[i0]+m[i3]+3*(m[i2]-m[i1]) );
 107       }
 108     }
 109   }
 110  }
 111
 112
 113
 114 inline Vc::float_v GetSpline3Vc(Vc::float_v v0, Vc::float_v v1, Vc::float_v v2, Vc::float_v v3,
 115                                 Vc::float_v x, Vc::float_v x2)
 116 {
 117   Vc::float_v dv = v2-v1;
 118   Vc::float_v z0 = 0.5f*(v2-v0);
 119   Vc::float_v z1 = 0.5f*(v3-v1);
 120   return x2*( (z1-dv + z0-dv)*(x-1) - (z0-dv) ) + z0*x + v1;
 121 }
 122
 123 void AliHLTTPCSpline2D3D::GetValue(Float_t A, Float_t B, Float_t XYZ[] ) const
 124 {
 125   //
 126   //  Get Interpolated value at A,B
 127   //
 128
 129   Float_t lA = (A-fMinA)*fScaleA-1.f;
 130   Int_t iA = (int) lA;
 131   if( lA<0 ) iA=0;
 132   else if( iA>fNA-4 ) iA = fNA-4;
 133
 134   Float_t lB = (B-fMinB)*fScaleB -1.f;
 135   Int_t iB = (int) lB;
 136   if( lB<0 ) iB=0;
 137   else if( iB>fNB-4 ) iB = fNB-4;
 138
 139   if( Vc::float_v::Size==4 ){
 140     Vc::float_v da = lA-iA;
 141     Vc::float_v db = lB-iB;
 142     Vc::float_v db2=db*db;
 143
 144     Vc::float_v v[4];
 145     Int_t ind = iA*fNB  + iB;
 146     const Vc::float_v *m = reinterpret_cast< const Vc::float_v *>(fXYZ);
 147
 148     for( Int_t i=0; i<4; i++ ){
 149       v[i] = GetSpline3Vc(m[ind+0],m[ind+1],m[ind+2],m[ind+3],db,db2);
 150       ind+=fNB;
 151     }
 152     Vc::float_v res=GetSpline3Vc(v[0],v[1],v[2],v[3],da,da*da);
 153     XYZ[0] =  res[0];
 154     XYZ[1] =  res[1];
 155     XYZ[2] =  res[2];
 156   } else {
 157     Float_t da = lA-iA;
 158     Float_t db = lB-iB;
 159
 160     Float_t vx[4];
 161     Float_t vy[4];
 162     Float_t vz[4];
 163     Int_t ind = iA*fNB  + iB;
 164     for( Int_t i=0; i<4; i++ ){
 165       int ind4 = ind*4;
 166       vx[i] = GetSpline3(fXYZ[ind4+0],fXYZ[ind4+4],fXYZ[ind4 +8],fXYZ[ind4+12],db);
 167       vy[i] = GetSpline3(fXYZ[ind4+1],fXYZ[ind4+5],fXYZ[ind4 +9],fXYZ[ind4+13],db);
 168       vz[i] = GetSpline3(fXYZ[ind4+2],fXYZ[ind4+6],fXYZ[ind4+10],fXYZ[ind4+14],db);
 169       ind+=fNB;
 170     }
 171     XYZ[0] =  GetSpline3(vx,da);
 172     XYZ[1] =  GetSpline3(vy,da);
 173     XYZ[2] =  GetSpline3(vz,da);
 174   }
 175 }