Line data Source code
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 : #include "AliHLTTPCSpline2D3D.h"
24 : #include "AliHLTTPCSpline2D3DObject.h"
25 : #include "Vc/Vc"
26 :
27 : #include <iostream>
28 : #include <iomanip>
29 :
30 : using namespace std;
31 :
32 : void AliHLTTPCSpline2D3D::Init(Float_t minA,Float_t maxA, Int_t nBinsA, Float_t minB,Float_t maxB, Int_t nBinsB)
33 : {
34 : //
35 : // Initialisation
36 : //
37 :
38 0 : if( maxA<= minA ) maxA = minA+1;
39 0 : if( maxB<= minB ) maxB = minB+1;
40 0 : if( nBinsA <4 ) nBinsA = 4;
41 0 : if( nBinsB <4 ) nBinsB = 4;
42 :
43 0 : fNA = nBinsA;
44 0 : fNB = nBinsB;
45 0 : fN = fNA*fNB;
46 0 : fMinA = minA;
47 0 : fMinB = minB;
48 :
49 0 : fStepA = (maxA-minA)/(nBinsA-1);
50 0 : fStepB = (maxB-minB)/(nBinsB-1);
51 0 : fScaleA = 1./fStepA;
52 0 : fScaleB = 1./fStepB;
53 :
54 0 : Vc::free( fXYZ );
55 0 : fXYZ = Vc::malloc< float, Vc::AlignOnCacheline>( 4*fN );
56 0 : memset ( fXYZ, 0, fN*4*sizeof(Float_t) );
57 0 : }
58 :
59 : AliHLTTPCSpline2D3D::~AliHLTTPCSpline2D3D()
60 0 : {
61 0 : Vc::free( fXYZ );
62 0 : }
63 :
64 : Int_t AliHLTTPCSpline2D3D::WriteToBuffer( char *buf, size_t &size ) const
65 : {
66 0 : size_t maxSize = size;
67 0 : size = 0;
68 0 : if( EstimateBufferSize() > maxSize ) return -1;
69 0 : Int_t *pInt = reinterpret_cast<Int_t*> (buf);
70 0 : pInt[0] = Version();
71 0 : pInt[1] = fNA;
72 0 : pInt[2] = fNB;
73 0 : Float_t *pFloat = reinterpret_cast<Float_t*> (pInt+3);
74 0 : pFloat[0] = fMinA;
75 0 : pFloat[1] = fMinB;
76 0 : pFloat[2] = fStepA;
77 0 : pFloat[3] = fStepB;
78 0 : pFloat+=4;
79 0 : for( int i=0; i<fN; i+=4, pFloat+=3 ){
80 0 : pFloat[0] = fXYZ[i+0];
81 0 : pFloat[1] = fXYZ[i+1];
82 0 : pFloat[2] = fXYZ[i+2];
83 : }
84 0 : size = EstimateBufferSize();
85 : return 0;
86 0 : }
87 :
88 : Int_t AliHLTTPCSpline2D3D::ReadFromBuffer( const char*buf, size_t &size )
89 : {
90 0 : size_t maxSize = size;
91 0 : size = 0;
92 0 : if( maxSize < sizeof(Int_t)*3 + sizeof(Float_t)*(4) ) return -1;
93 0 : const Int_t *pInt = reinterpret_cast< const Int_t*> (buf);
94 0 : const Float_t *pFloat = reinterpret_cast< const Float_t*> (pInt+3);
95 0 : if( pInt[0]!= Version() ) return -1;
96 0 : Int_t nA = pInt[1];
97 0 : Int_t nB = pInt[2];
98 0 : Int_t n = nA*nB;
99 0 : if( maxSize < sizeof(Int_t)*3 + sizeof(Float_t)*(4+3*n) ) return -2;
100 :
101 0 : fNA = nA;
102 0 : fNB = nB;
103 0 : fN = n;
104 0 : fMinA = pFloat[0];
105 0 : fMinB = pFloat[1];
106 0 : fStepA = pFloat[2];
107 0 : fStepB = pFloat[3];
108 : int err = 0;
109 0 : if( fabs( fStepA ) <1.e-8 ) { fStepA = 1.; err=-3; }
110 0 : if( fabs( fStepB ) <1.e-8 ) { fStepB = 1.; err=-3; }
111 :
112 0 : fScaleA = 1./fStepA;
113 0 : fScaleB = 1./fStepB;
114 :
115 0 : pFloat+=4;
116 0 : Vc::free( fXYZ );
117 0 : fXYZ = Vc::malloc< float, Vc::AlignOnCacheline>( 4*fN );
118 0 : for( int i=0; i<4*fN; i+=4, pFloat+=3 ){
119 0 : fXYZ[i+0] = pFloat[0];
120 0 : fXYZ[i+1] = pFloat[1];
121 0 : fXYZ[i+2] = pFloat[2];
122 0 : fXYZ[i+3] = 0;
123 : }
124 0 : size = sizeof(Int_t)*3 + sizeof(Float_t)*(4+3*n);
125 : return err;
126 0 : }
127 :
128 : void AliHLTTPCSpline2D3D::WriteToObject( AliHLTTPCSpline2D3DObject &obj )
129 : {
130 : //
131 : // write spline to ROOT object to store it in database
132 : //
133 0 : obj.SetMinA( fMinA );
134 0 : obj.SetMinB( fMinB );
135 0 : obj.SetStepA( fStepA );
136 0 : obj.SetStepB( fStepB );
137 0 : obj.InitGrid( fNA, fNB );
138 0 : for( int i=0; i<fN; i++ ) obj.SetGridValue( i, fXYZ + i*4);
139 0 : }
140 :
141 : void AliHLTTPCSpline2D3D::ReadFromObject( const AliHLTTPCSpline2D3DObject &obj )
142 : {
143 : //
144 : // read spline from ROOT object stored in database
145 : //
146 :
147 0 : fMinA = obj.GetMinA();
148 0 : fStepA = obj.GetStepA();
149 0 : fNA = obj.GetNPointsA();
150 :
151 0 : fMinB = obj.GetMinB();
152 0 : fStepB = obj.GetStepB();
153 0 : fNB = obj.GetNPointsB();
154 :
155 0 : fN = fNA*fNB;
156 :
157 0 : fScaleA = fabs(fStepA) >1.e-8 ?1./fStepA :1.;
158 0 : fScaleB = fabs(fStepB) >1.e-8 ?1./fStepB :1.;
159 :
160 0 : Vc::free( fXYZ );
161 0 : fXYZ = Vc::malloc< float, Vc::AlignOnCacheline>( 4*fN );
162 0 : for( int i=0; i<fN; i++ ){
163 0 : obj.GetGridValue( i, fXYZ + i*4 );
164 0 : fXYZ[i*4+3] = 0.;
165 : }
166 0 : }
167 :
168 : void AliHLTTPCSpline2D3D::Consolidate()
169 : {
170 : //
171 : // Consolidate the map
172 : //
173 :
174 0 : Float_t *m = fXYZ;
175 0 : for( int iXYZ=0; iXYZ<3; iXYZ++ ){
176 0 : for( int iA=0; iA<fNA; iA++){
177 : {
178 0 : int i0 = 4*iA*fNB + iXYZ;
179 0 : int i1 = i0+4;
180 0 : int i2 = i0+4*2;
181 0 : int i3 = i0+4*3;
182 0 : m[i0] = 0.5*( m[i3]+m[i0]+3*(m[i1]-m[i2]) );
183 : }
184 : {
185 0 : int i0 = 4*(iA*fNB+fNB-4) + iXYZ;
186 0 : int i1 = i0+4;
187 0 : int i2 = i0+4*2;
188 0 : int i3 = i0+4*3;
189 0 : m[i3] = 0.5*( m[i0]+m[i3]+3*(m[i2]-m[i1]) );
190 : }
191 : }
192 0 : for( int iB=0; iB<fNB; iB++){
193 : {
194 0 : int i0 = 4*(0*fNB +iB) +iXYZ;
195 0 : int i1 = i0+4*fNB;
196 0 : int i2 = i1+4*fNB;
197 0 : int i3 = i2+4*fNB;
198 0 : m[i0] = 0.5*( m[i3]+m[i0]+3*(m[i1]-m[i2]) );
199 : }
200 : {
201 0 : int i0 = 4*((fNA-4)*fNB +iB) +iXYZ;
202 0 : int i1 = i0+4*fNB;
203 0 : int i2 = i1+4*fNB;
204 0 : int i3 = i2+4*fNB;
205 0 : m[i3] = 0.5*( m[i0]+m[i3]+3*(m[i2]-m[i1]) );
206 : }
207 : }
208 : }
209 0 : }
210 :
211 :
212 :
213 : inline Vc::float_v GetSpline3Vc(Vc::float_v v0, Vc::float_v v1, Vc::float_v v2, Vc::float_v v3,
214 : Vc::float_v x, Vc::float_v x2)
215 : {
216 0 : Vc::float_v dv = v2-v1;
217 0 : Vc::float_v z0 = 0.5f*(v2-v0);
218 0 : Vc::float_v z1 = 0.5f*(v3-v1);
219 0 : return x2*( (z1-dv + z0-dv)*(x-1) - (z0-dv) ) + z0*x + v1;
220 0 : }
221 :
222 : void AliHLTTPCSpline2D3D::GetValue(Float_t A, Float_t B, Float_t XYZ[] ) const
223 : {
224 : //
225 : // Get Interpolated value at A,B
226 : //
227 :
228 0 : Float_t lA = (A-fMinA)*fScaleA-1.f;
229 0 : Int_t iA = (int) lA;
230 0 : if( lA<0 ) iA=0;
231 0 : else if( iA>fNA-4 ) iA = fNA-4;
232 :
233 0 : Float_t lB = (B-fMinB)*fScaleB -1.f;
234 0 : Int_t iB = (int) lB;
235 0 : if( lB<0 ) iB=0;
236 0 : else if( iB>fNB-4 ) iB = fNB-4;
237 :
238 : if( Vc::float_v::Size==4 ){
239 0 : Vc::float_v da = lA-iA;
240 0 : Vc::float_v db = lB-iB;
241 0 : Vc::float_v db2=db*db;
242 :
243 0 : Vc::float_v v[4];
244 0 : Int_t ind = iA*fNB + iB;
245 0 : const Vc::float_v *m = reinterpret_cast< const Vc::float_v *>(fXYZ);
246 :
247 0 : for( Int_t i=0; i<4; i++ ){
248 0 : v[i] = GetSpline3Vc(m[ind+0],m[ind+1],m[ind+2],m[ind+3],db,db2);
249 0 : ind+=fNB;
250 : }
251 0 : Vc::float_v res=GetSpline3Vc(v[0],v[1],v[2],v[3],da,da*da);
252 0 : XYZ[0] = res[0];
253 0 : XYZ[1] = res[1];
254 0 : XYZ[2] = res[2];
255 0 : } else {
256 : Float_t da = lA-iA;
257 : Float_t db = lB-iB;
258 :
259 : Float_t vx[4];
260 : Float_t vy[4];
261 : Float_t vz[4];
262 : Int_t ind = iA*fNB + iB;
263 : for( Int_t i=0; i<4; i++ ){
264 : int ind4 = ind*4;
265 : vx[i] = GetSpline3(fXYZ[ind4+0],fXYZ[ind4+4],fXYZ[ind4 +8],fXYZ[ind4+12],db);
266 : vy[i] = GetSpline3(fXYZ[ind4+1],fXYZ[ind4+5],fXYZ[ind4 +9],fXYZ[ind4+13],db);
267 : vz[i] = GetSpline3(fXYZ[ind4+2],fXYZ[ind4+6],fXYZ[ind4+10],fXYZ[ind4+14],db);
268 : ind+=fNB;
269 : }
270 : XYZ[0] = GetSpline3(vx,da);
271 : XYZ[1] = GetSpline3(vy,da);
272 : XYZ[2] = GetSpline3(vz,da);
273 : }
274 0 : }
|
