Line data Source code
1 : // $Id$
2 :
3 : //**************************************************************************
4 : //* This file is property of and copyright by the ALICE HLT Project *
5 : //* ALICE Experiment at CERN, All rights reserved. *
6 : //* *
7 : //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
8 : //* for The ALICE HLT Project. *
9 : //* *
10 : //* Permission to use, copy, modify and distribute this software and its *
11 : //* documentation strictly for non-commercial purposes is hereby granted *
12 : //* without fee, provided that the above copyright notice appears in all *
13 : //* copies and that both the copyright notice and this permission notice *
14 : //* appear in the supporting documentation. The authors make no claims *
15 : //* about the suitability of this software for any purpose. It is *
16 : //* provided "as is" without express or implied warranty. *
17 : //**************************************************************************
18 :
19 : /// @file AliHLTTPCTrackGeometry.cxx
20 : /// @author Matthias Richter
21 : /// @date 2011-05-20
22 : /// @brief Desciption of a track by a sequence of track points
23 : ///
24 :
25 : #include "AliHLTTPCTrackGeometry.h"
26 : #include "AliHLTTPCGeometry.h"
27 : #include "AliHLTTPCSpacePointData.h"
28 : #include "AliHLTTPCClusterDataFormat.h"
29 : #include "AliHLTTPCSpacePointContainer.h"
30 : #include "AliHLTTPCRawSpacePointContainer.h"
31 : #include "AliHLTTPCDefinitions.h"
32 : #include "AliHLTComponent.h"
33 : #include "AliHLTGlobalBarrelTrack.h"
34 : #include "AliHLTDataDeflater.h"
35 : #include "AliHLTErrorGuard.h"
36 : #include "TMath.h"
37 : #include "TH2F.h"
38 : #include <memory>
39 : #include <sstream>
40 : using namespace std;
41 :
42 : /** ROOT macro for the implementation of ROOT specific class methods */
43 6 : ClassImp(AliHLTTPCTrackGeometry)
44 :
45 : AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry()
46 0 : : AliHLTTrackGeometry()
47 0 : , fRawTrackPoints()
48 0 : , fDriftTimeFactorA(0.)
49 0 : , fDriftTimeOffsetA(0.)
50 0 : , fDriftTimeFactorC(0.)
51 0 : , fDriftTimeOffsetC(0.)
52 0 : {
53 : /// standard constructor
54 0 : }
55 :
56 : AliHLTTPCTrackGeometry::AliHLTTPCTrackGeometry(const AliHLTTPCTrackGeometry& src)
57 0 : : AliHLTTrackGeometry(src)
58 0 : , fRawTrackPoints(src.fRawTrackPoints)
59 0 : , fDriftTimeFactorA(0.)
60 0 : , fDriftTimeOffsetA(0.)
61 0 : , fDriftTimeFactorC(0.)
62 0 : , fDriftTimeOffsetC(0.)
63 0 : {
64 : /// copy constructor
65 0 : }
66 :
67 : AliHLTTPCTrackGeometry& AliHLTTPCTrackGeometry::operator=(const AliHLTTPCTrackGeometry& src)
68 : {
69 : /// assignment operator
70 0 : AliHLTTrackGeometry::operator=(src);
71 0 : fRawTrackPoints.assign(src.fRawTrackPoints.begin(), src.fRawTrackPoints.end());
72 0 : return *this;
73 : }
74 :
75 0 : AliHLTTPCTrackGeometry::~AliHLTTPCTrackGeometry()
76 0 : {
77 : /// destructor
78 0 : }
79 :
80 : float AliHLTTPCTrackGeometry::GetPlaneAlpha(AliHLTUInt32_t planeId) const
81 : {
82 : /// alpha of the plane
83 0 : UInt_t slice=AliHLTTPCSpacePointData::GetSlice(planeId);
84 0 : float alpha=( slice + 0.5 ) * TMath::Pi() / 9.0;
85 0 : if (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
86 0 : return alpha;
87 : }
88 :
89 : float AliHLTTPCTrackGeometry::GetPlaneR(AliHLTUInt32_t planeId) const
90 : {
91 : /// radial distance from global {0,0,0}
92 0 : UInt_t partition=AliHLTTPCSpacePointData::GetPatch(planeId);
93 0 : UInt_t number=AliHLTTPCSpacePointData::GetNumber(planeId);
94 0 : Int_t row=AliHLTTPCGeometry::GetFirstRow(partition)+number;
95 0 : return AliHLTTPCGeometry::Row2X(row);
96 : }
97 :
98 : float AliHLTTPCTrackGeometry::GetPlaneTheta(AliHLTUInt32_t /*planeId*/) const
99 : {
100 : /// theta of the plane
101 0 : return 0.0;
102 : }
103 :
104 : bool AliHLTTPCTrackGeometry::CheckBounds(AliHLTUInt32_t planeId, float u, float /*v*/) const
105 : {
106 : /// check bounds in u and v coordinate
107 0 : float r=GetPlaneR(planeId);
108 0 : if (r<AliHLTTPCGeometry::GetFirstRow(0)) return false;
109 :
110 : // TODO: check if the pad width needs to be considered here
111 0 : return TMath::Abs(TMath::ASin(u/r))<=TMath::Pi()/18;
112 0 : }
113 :
114 : int AliHLTTPCTrackGeometry::CalculateTrackPoints(const AliHLTExternalTrackParam& track)
115 : {
116 : /// calculate the track points, expects the global magnetic field to be initialized
117 0 : AliHLTGlobalBarrelTrack bt(track);
118 0 : return CalculateTrackPoints(bt);
119 0 : }
120 :
121 : int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track)
122 : {
123 : /// calculate the track points, expects the global magnetic field to be initialized
124 : /// depending on the x coordinate of the first track point the corresponding padrow
125 : /// is searched and the points are calculated outwards. Eventually the points are
126 : /// calculated inwards as a second step.
127 : int iResult=0;
128 : int firstpadrow=0;
129 0 : for (;
130 0 : firstpadrow<AliHLTTPCGeometry::GetNRows() &&
131 0 : AliHLTTPCGeometry::Row2X(firstpadrow)+AliHLTTPCGeometry::GetPadLength(firstpadrow)<track.GetX();
132 0 : firstpadrow++);
133 0 : if (firstpadrow>=AliHLTTPCGeometry::GetNRows()) return 0;
134 : // first calculated outwards
135 0 : iResult=CalculateTrackPoints(track, firstpadrow, 1);
136 0 : if (iResult>=0 && firstpadrow>0)
137 : // now calculate inwards
138 0 : iResult=CalculateTrackPoints(track, firstpadrow-1, -1);
139 0 : return iResult;
140 0 : }
141 :
142 : int AliHLTTPCTrackGeometry::CalculateTrackPoints(AliHLTGlobalBarrelTrack& track, int firstpadrow, int step)
143 : {
144 : /// calculate the track points, expects the global magnetic field to be initialized
145 : float offsetAlpha=0.0;
146 0 : for (int padrow=firstpadrow; padrow>=0 && padrow<AliHLTTPCGeometry::GetNRows(); padrow+=step) {
147 0 : float x=AliHLTTPCGeometry::Row2X(padrow);
148 0 : float y=0.0;
149 0 : float z=0.0;
150 :
151 : int maxshift=9;
152 : int shift=0;
153 : int result=0;
154 0 : do {
155 : // start calculation of crossing points with padrow planes in the slice of the first point
156 : // plane alpha corresponds to alpha of the track, switch to neighboring slice if the result
157 : // is out of bounds
158 0 : if ((result=track.CalculateCrossingPoint(x, track.GetAlpha()-offsetAlpha, y, z))<1) break;
159 0 : float pointAlpha=TMath::ATan(y/x);
160 0 : if (TMath::Abs(pointAlpha)>TMath::Pi()/18) {
161 0 : offsetAlpha+=(pointAlpha>0?-1:1)*TMath::Pi()/9;
162 : result=0;
163 0 : }
164 0 : } while (result==0 && shift++<maxshift);
165 0 : if (result<1) continue;
166 0 : float planealpha=track.GetAlpha()-offsetAlpha;
167 0 : if (planealpha<0) planealpha+=TMath::TwoPi();
168 0 : if (planealpha>TMath::TwoPi()) planealpha-=TMath::TwoPi();
169 0 : int slice=int(9*planealpha/TMath::Pi());
170 0 : if (z<0) slice+=18;
171 0 : if (slice>=36) {
172 0 : HLTError("invalid slice %d calculated from alpha %f", slice, track.GetAlpha());
173 : }
174 0 : int partition=AliHLTTPCGeometry::GetPatch(padrow);
175 0 : int row=padrow-AliHLTTPCGeometry::GetFirstRow(partition);
176 0 : UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
177 0 : if (TMath::Abs(planealpha-GetPlaneAlpha(id))>0.0001) {
178 0 : HLTError("alpha missmatch for plane %08x (slice %d): alpha from id %f (%.0f deg), expected %f (%.0f deg)", id, slice, GetPlaneAlpha(id), 180*GetPlaneAlpha(id)/TMath::Pi(), planealpha, 180*planealpha/TMath::Pi());
179 : }
180 0 : if (AddTrackPoint(AliHLTTrackPoint(id, y, z), AliHLTTPCSpacePointData::GetID(slice, partition, 0))>=0) {
181 0 : Float_t rpt[3]={0.,y,z}; // row pad time
182 0 : AliHLTTPCGeometry::LocHLT2Raw(rpt, slice, padrow);
183 0 : float m=fDriftTimeFactorA;
184 0 : float n=fDriftTimeOffsetA;
185 0 : if (slice>=18) {
186 0 : m=fDriftTimeFactorC;
187 0 : n=fDriftTimeOffsetC;
188 0 : }
189 0 : if (TMath::Abs(m)>0.) {
190 0 : rpt[2]=(z-n)/m;
191 0 : if (step>0) {
192 : // ascending padrows added at end
193 0 : fRawTrackPoints.push_back(AliHLTTrackPoint(id, rpt[1], rpt[2]));
194 0 : } else {
195 : // descending padrows added at begin
196 0 : fRawTrackPoints.insert(fRawTrackPoints.begin(), AliHLTTrackPoint(id, rpt[1], rpt[2]));
197 : }
198 : // FIXME: implement a Print function for the raw track points
199 : // stringstream sout;
200 : // sout << " slice " << setfill(' ') << setw(3) << fixed << right << slice
201 : // << " row " << setfill(' ') << setw(3) << fixed << right << padrow
202 : // << " id 0x" << hex << setw(8) << id << dec
203 : // << " y " << setfill(' ') << setw(5) << fixed << right << setprecision (2) << y
204 : // << " z " << setfill(' ') << setw(5) << fixed << right << setprecision (2) << z
205 : // << " pad " << setfill(' ') << setw(5) << fixed << right << setprecision (2) << rpt[1]
206 : // << " time " << setfill(' ') << setw(5) << fixed << right << setprecision (2) << rpt[2]
207 : // << endl;
208 : // cout << sout.str();
209 : } else {
210 0 : ALIHLTERRORGUARD(1, "drift time correction not initialized, can not add track points in raw coordinates");
211 : }
212 0 : }
213 0 : }
214 0 : return 0;
215 0 : }
216 :
217 : int AliHLTTPCTrackGeometry::FindMatchingTrackPoint(AliHLTUInt32_t spacepointId, float spacepoint[3], AliHLTUInt32_t& planeId)
218 : {
219 : /// find the track point which can be associated to a spacepoint with coordinates and id
220 0 : UInt_t slice=AliHLTTPCSpacePointData::GetSlice(spacepointId);
221 0 : UInt_t partition=AliHLTTPCSpacePointData::GetPatch(spacepointId);
222 0 : int row=AliHLTTPCGeometry::GetPadRow(spacepoint[0]);
223 0 : bool bSpecialRow=row==30 || row==90 || row==139;
224 0 : if (row<AliHLTTPCGeometry::GetFirstRow(partition) || row>AliHLTTPCGeometry::GetLastRow(partition)) {
225 0 : HLTError("row number %d calculated from x value %f is outside slice %d partition %d", row, spacepoint[0], slice, partition);
226 0 : return -EINVAL;
227 : }
228 :
229 : // find the crossing point of the track with the padrow plane where
230 : // the spacepoint is
231 : // 1) calculate plane id from slice, partition and row (within partition)
232 0 : row-=AliHLTTPCGeometry::GetFirstRow(partition);
233 0 : UInt_t id=AliHLTTPCSpacePointData::GetID(slice, partition, row);
234 0 : const AliHLTTrackPoint* point=GetTrackPoint(id);
235 : // track might be outside the partition and cross the central membrane
236 : // search in the other half of the TPC
237 0 : if (!point && slice<18) {
238 : // search in the neighboring partition on the C side
239 0 : id=AliHLTTPCSpacePointData::GetID(slice+18, partition, row);
240 0 : point=GetTrackPoint(id);
241 0 : } else if (!point && slice>=18) {
242 : // search in the neighboring partition on the A side
243 0 : id=AliHLTTPCSpacePointData::GetID(slice-18, partition, row);
244 0 : point=GetTrackPoint(id);
245 0 : }
246 :
247 : // search in the neighboring partition, this takes account for rows
248 : // 30, 90, and 139 which are partly in one and the other partition
249 0 : if (!point && bSpecialRow) {
250 0 : row+=AliHLTTPCGeometry::GetFirstRow(partition);
251 0 : row-=AliHLTTPCGeometry::GetFirstRow(partition-1);
252 0 : id=AliHLTTPCSpacePointData::GetID(slice, partition-1, row);
253 0 : point=GetTrackPoint(id);
254 0 : if (!point && slice<18) {
255 : // search in the neighboring partition on the C side
256 0 : id=AliHLTTPCSpacePointData::GetID(slice+18, partition-1, row);
257 0 : point=GetTrackPoint(id);
258 0 : } else if (!point && slice>=18) {
259 : // search in the neighboring partition on the A side
260 0 : id=AliHLTTPCSpacePointData::GetID(slice-18, partition-1, row);
261 0 : point=GetTrackPoint(id);
262 0 : }
263 : }
264 :
265 0 : if (point) {
266 0 : planeId=id;
267 0 : if (point->HaveAssociatedSpacePoint()) {
268 0 : if (GetVerbosity()>2) HLTInfo("descarding spacepoint 0x%08x z=%f y=%f z=%f: track point 0x%08x already occupied", spacepoint[0], spacepoint[1], spacepoint[2], planeId);
269 0 : return 0; // already occupied
270 : }
271 : float maxdy=2.;
272 : float maxdz=2.;
273 0 : if (TMath::Abs(point->GetU()-spacepoint[1])>maxdy) {
274 0 : if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x y %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetU(), maxdy);
275 0 : return -ENOENT;
276 : }
277 0 : if (TMath::Abs(point->GetV()-spacepoint[2])>maxdz) {
278 0 : if (GetVerbosity()>0) HLTInfo("descarding spacepoint 0x%08x y=%f z=%f: track point 0x%08x z %f outside tolerance %f", spacepoint[1], spacepoint[2], planeId, point->GetV(), maxdz);
279 0 : return -ENOENT;
280 : }
281 0 : return 1;
282 : }
283 0 : return -ENOENT;
284 0 : }
285 :
286 :
287 : int AliHLTTPCTrackGeometry::RegisterTrackPoints(AliHLTTrackGrid* pGrid) const
288 : {
289 : /// register track points in the index grid, at this step the number
290 : /// of tracks in each cell is counted
291 0 : if (!pGrid) return -EINVAL;
292 : int iResult=0;
293 0 : for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
294 0 : tp!=TrackPoints().end() && iResult>=0; tp++) {
295 0 : AliHLTUInt32_t id=tp->GetId();
296 0 : iResult=pGrid->CountSpacePoint(AliHLTTPCSpacePointData::GetSlice(id),
297 0 : AliHLTTPCSpacePointData::GetPatch(id),
298 0 : AliHLTTPCSpacePointData::GetNumber(id));
299 : }
300 : return iResult;
301 0 : }
302 :
303 : int AliHLTTPCTrackGeometry::FillTrackPoints(AliHLTTrackGrid* pGrid) const
304 : {
305 : /// fill track points to index grid
306 0 : if (!pGrid) return -EINVAL;
307 : int iResult=0;
308 0 : for (vector<AliHLTTrackPoint>::const_iterator tp=TrackPoints().begin();
309 0 : tp!=TrackPoints().end() && iResult>=0; tp++) {
310 0 : AliHLTUInt32_t id=tp->GetId();
311 0 : iResult=pGrid->AddSpacePoint(GetTrackId(),
312 0 : AliHLTTPCSpacePointData::GetSlice(id),
313 0 : AliHLTTPCSpacePointData::GetPatch(id),
314 0 : AliHLTTPCSpacePointData::GetNumber(id));
315 : }
316 : return iResult;
317 0 : }
318 :
319 : AliHLTSpacePointContainer* AliHLTTPCTrackGeometry::ConvertToSpacePoints(bool bAssociated) const
320 : {
321 : /// create a collection of all points
322 0 : std::auto_ptr<AliHLTTPCSpacePointContainer> spacepoints(new AliHLTTPCSpacePointContainer);
323 0 : if (!spacepoints.get()) return NULL;
324 :
325 0 : const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
326 : unsigned i=0;
327 0 : while (i<trackPoints.size()) {
328 : // allocate buffer for all points, even though the buffer might not be filled
329 : // completely because of a partition change
330 0 : int nofPoints=trackPoints.size()-i;
331 0 : int blocksize=sizeof(AliHLTTPCClusterData)+nofPoints*sizeof(AliHLTTPCSpacePointData);
332 0 : AliHLTUInt8_t* pBuffer=spacepoints->Alloc(blocksize);
333 0 : if (!pBuffer) return NULL;
334 0 : AliHLTTPCClusterData* pClusterData=reinterpret_cast<AliHLTTPCClusterData*>(pBuffer);
335 0 : pClusterData->fSpacePointCnt=0;
336 0 : AliHLTTPCSpacePointData* pClusters=pClusterData->fSpacePoints;
337 : int currentSlice=-1;
338 : int currentPartition=-1;
339 0 : for (; i<trackPoints.size(); i++) {
340 0 : if (bAssociated && !trackPoints[i].HaveAssociatedSpacePoint()) continue;
341 0 : AliHLTUInt32_t planeId=trackPoints[i].GetId();
342 0 : int slice=AliHLTTPCSpacePointData::GetSlice(planeId);
343 0 : int partition=AliHLTTPCSpacePointData::GetPatch(planeId);
344 0 : int number=AliHLTTPCSpacePointData::GetNumber(planeId);
345 : if ((currentSlice>=0 && currentSlice!=slice) || (currentPartition>=0 && currentPartition!=partition)) {
346 : // change of partition or slice, need to go to next block
347 : // 2011-07-26 currently all spacepoints go into one block, if separated
348 : // blocks per partition are needed one has to leave the inner loop here
349 : // and set the data block specification below
350 : // Caution: not tested, only the last block seems to make it through
351 : //break;
352 : }
353 : currentSlice=slice;
354 : currentPartition=partition;
355 0 : pClusters[pClusterData->fSpacePointCnt].fX=GetPlaneR(planeId);
356 0 : pClusters[pClusterData->fSpacePointCnt].fY=trackPoints[i].GetU();
357 0 : pClusters[pClusterData->fSpacePointCnt].fZ=trackPoints[i].GetV();
358 0 : pClusters[pClusterData->fSpacePointCnt].fID=planeId;
359 0 : pClusters[pClusterData->fSpacePointCnt].fPadRow=AliHLTTPCGeometry::GetFirstRow(partition)+number;
360 0 : pClusters[pClusterData->fSpacePointCnt].fSigmaY2=0.;
361 0 : pClusters[pClusterData->fSpacePointCnt].fSigmaZ2=0.;
362 0 : pClusters[pClusterData->fSpacePointCnt].fCharge=0;
363 0 : pClusters[pClusterData->fSpacePointCnt].fQMax=0;
364 0 : pClusters[pClusterData->fSpacePointCnt].fUsed=0;
365 0 : pClusters[pClusterData->fSpacePointCnt].fTrackN=0;
366 0 : pClusterData->fSpacePointCnt++;
367 0 : }
368 0 : AliHLTComponentBlockData bd;
369 0 : AliHLTComponent::FillBlockData(bd);
370 0 : bd.fPtr=pBuffer;
371 0 : bd.fSize=sizeof(AliHLTTPCClusterData)+pClusterData->fSpacePointCnt*sizeof(AliHLTTPCSpacePointData);
372 0 : AliHLTComponent::SetDataType(bd.fDataType, "CLUSTERS", "TPC ");
373 0 : bd.fSpecification=kAliHLTVoidDataSpec;//AliHLTTPCDefinitions::EncodeDataSpecification(currentSlice, currentSlice, currentPartition, currentPartition);
374 0 : spacepoints->AddInputBlock(&bd);
375 0 : }
376 :
377 0 : return spacepoints.release();
378 0 : }
379 :
380 : const AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id) const
381 : {
382 : /// get raw track point of id
383 0 : const AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
384 0 : if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
385 0 : return p;
386 0 : }
387 :
388 : AliHLTTrackGeometry::AliHLTTrackPoint* AliHLTTPCTrackGeometry::GetRawTrackPoint(AliHLTUInt32_t id)
389 : {
390 : /// get raw track point of id
391 0 : AliHLTTrackGeometry::AliHLTTrackPoint* p=find(&fRawTrackPoints[0], &fRawTrackPoints[fRawTrackPoints.size()], id);
392 0 : if (p==&fRawTrackPoints[fRawTrackPoints.size()]) return 0;
393 0 : return p;
394 0 : }
395 :
396 : int AliHLTTPCTrackGeometry::FillRawResidual(int coordinate, TH2* histo, AliHLTSpacePointContainer* points) const
397 : {
398 : // fill residual histogram
399 0 : if (!histo || !points) return -EINVAL;
400 0 : const vector<AliHLTTrackPoint>& trackPoints=TrackPoints();
401 0 : for (vector<AliHLTTrackPoint>::const_iterator trackpoint=trackPoints.begin();
402 0 : trackpoint!=trackPoints.end(); trackpoint++) {
403 0 : if (!trackpoint->HaveAssociatedSpacePoint()) continue;
404 0 : for (vector<AliHLTTrackSpacepoint>::const_iterator sp=(trackpoint->GetSpacepoints()).begin();
405 0 : sp!=(trackpoint->GetSpacepoints()).end(); sp++) {
406 0 : AliHLTUInt32_t spacepointId=sp->fId;
407 0 : vector<AliHLTTrackPoint>::const_iterator rawpoint=find(fRawTrackPoints.begin(), fRawTrackPoints.end(), trackpoint->GetId());
408 0 : if (rawpoint==fRawTrackPoints.end()) {
409 0 : HLTError("can not find track raw coordinates of track point 0x%08x", trackpoint->GetId());
410 0 : continue;
411 : }
412 0 : if (!points->Check(spacepointId)) {
413 : //HLTError("can not find associated space point 0x%08x of track point 0x%08x", spacepointId, trackpoint->GetId());
414 0 : continue;
415 : }
416 : float value=0.;
417 0 : if (coordinate==0) {
418 0 : value=rawpoint->GetU()-points->GetY(spacepointId);
419 0 : histo->Fill(GetPlaneR(trackpoint->GetId()), value);
420 0 : } else {
421 0 : value=rawpoint->GetV()-points->GetZ(spacepointId);
422 : //histo->Fill(GetPlaneR(trackpoint->GetId()), value);
423 0 : histo->Fill(rawpoint->GetV(), value);
424 : }
425 0 : }
426 0 : }
427 : return 0;
428 0 : }
429 :
430 : int AliHLTTPCTrackGeometry::Write(const AliHLTGlobalBarrelTrack& track,
431 : AliHLTSpacePointContainer* pSpacePoints,
432 : AliHLTDataDeflater* pDeflater,
433 : AliHLTUInt8_t* outputPtr,
434 : AliHLTUInt32_t size,
435 : vector<AliHLTUInt32_t>* writtenClusterIds,
436 : const char* option) const
437 : {
438 : // write track block to buffer
439 0 : if (size<=sizeof(AliHLTTPCTrackBlock)) return -ENOSPC;
440 0 : AliHLTTPCTrackBlock* pTrackBlock=reinterpret_cast<AliHLTTPCTrackBlock*>(outputPtr);
441 0 : pTrackBlock->fSize=sizeof(AliHLTTPCTrackBlock); // size of cluster block added later
442 0 : float alpha=track.GetAlpha();
443 0 : while (alpha<0.) alpha+=TMath::TwoPi();
444 0 : while (alpha>TMath::TwoPi()) alpha-=TMath::TwoPi();
445 0 : pTrackBlock->fSlice=AliHLTUInt8_t(9*alpha/TMath::Pi());
446 0 : if (pTrackBlock->fSlice>=36) {
447 0 : HLTError("invalid slice %d calculated from alpha %f", pTrackBlock->fSlice, track.GetAlpha());
448 : }
449 0 : pTrackBlock->fReserved=0;
450 0 : pTrackBlock->fX = track.GetX();
451 0 : pTrackBlock->fY = track.GetY();
452 0 : pTrackBlock->fZ = track.GetZ();
453 0 : pTrackBlock->fSinPsi = track.GetSnp();
454 0 : pTrackBlock->fTgl = track.GetTgl();
455 0 : pTrackBlock->fq1Pt = track.GetSigned1Pt();
456 :
457 0 : pDeflater->Clear();
458 0 : pDeflater->InitBitDataOutput(reinterpret_cast<AliHLTUInt8_t*>(outputPtr+sizeof(AliHLTTPCTrackBlock)), size-sizeof(AliHLTTPCTrackBlock));
459 0 : int result=WriteAssociatedClusters(pSpacePoints, pDeflater, writtenClusterIds, option);
460 0 : if (result<0) return result;
461 0 : pTrackBlock->fSize+=result;
462 0 : return pTrackBlock->fSize;
463 0 : }
464 :
465 : int AliHLTTPCTrackGeometry::WriteAssociatedClusters(AliHLTSpacePointContainer* pSpacePoints,
466 : AliHLTDataDeflater* pDeflater,
467 : vector<AliHLTUInt32_t>* writtenClusterIds,
468 : const char* /*option*/) const
469 : {
470 : // write associated clusters to buffer via deflater
471 0 : if (!pDeflater || !pSpacePoints) return -EINVAL;
472 0 : AliHLTTPCRawSpacePointContainer* pTPCRawSpacePoints=dynamic_cast<AliHLTTPCRawSpacePointContainer*>(pSpacePoints);
473 0 : if (!pTPCRawSpacePoints) {
474 0 : HLTError("invalid type of SpacePointContainer \"%s\", required AliHLTTPCRawSpacePointContainer", pSpacePoints->ClassName());
475 0 : return -EINVAL;
476 : }
477 : bool bReverse=true;
478 : bool bWriteSuccess=true;
479 : int writtenClusters=0;
480 : // filling of track points starts from first point on track outwards, and
481 : // then from that point inwards. That's why the lower padrows might be in
482 : // reverse order at the end of the track point array. If the last element
483 : // is bigger than the first element, only trackpoints in ascending order
484 : // are in the array
485 0 : vector<AliHLTTrackPoint>::const_iterator clrow=fRawTrackPoints.end();
486 0 : if (clrow!=fRawTrackPoints.begin()) {
487 0 : clrow--;
488 0 : AliHLTUInt32_t partition=AliHLTTPCSpacePointData::GetPatch(clrow->GetId());
489 0 : AliHLTUInt32_t partitionrow=AliHLTTPCSpacePointData::GetNumber(clrow->GetId());
490 0 : partitionrow+=AliHLTTPCGeometry::GetFirstRow(partition);
491 0 : AliHLTUInt32_t firstpartition=AliHLTTPCSpacePointData::GetPatch(fRawTrackPoints.begin()->GetId());
492 0 : AliHLTUInt32_t firstpartitionrow=AliHLTTPCSpacePointData::GetNumber(fRawTrackPoints.begin()->GetId());
493 0 : firstpartitionrow+=AliHLTTPCGeometry::GetFirstRow(firstpartition);
494 0 : if (partitionrow>=firstpartitionrow) {
495 : bReverse=false;
496 0 : clrow=fRawTrackPoints.begin();
497 0 : }
498 0 : }
499 0 : const AliHLTUInt32_t clusterCountBitLength=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kClusterCount].fBitLength;
500 0 : unsigned long dataPosition=pDeflater->GetCurrentByteOutputPosition();
501 0 : for (unsigned row=0; row<159 && bWriteSuccess; row++) {
502 0 : if (clrow!=fRawTrackPoints.end()) {
503 0 : AliHLTUInt32_t thisPartition=AliHLTTPCSpacePointData::GetPatch(clrow->GetId());
504 0 : AliHLTUInt32_t thisTrackRow=AliHLTTPCSpacePointData::GetNumber(clrow->GetId());
505 0 : thisTrackRow+=AliHLTTPCGeometry::GetFirstRow(thisPartition);
506 0 : if (thisTrackRow==row) {
507 : // write clusters
508 0 : const vector<AliHLTTrackSpacepoint>& clusters=clrow->GetSpacepoints();
509 0 : AliHLTUInt32_t haveClusters=clusters.size()>0;
510 : // 1 bit for clusters on that padrow
511 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(haveClusters);
512 0 : if (haveClusters) {
513 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputBits(clusters.size(), clusterCountBitLength);
514 0 : for (vector<AliHLTTrackSpacepoint>::const_iterator clid=clusters.begin();
515 0 : clid!=clusters.end() && bWriteSuccess; clid++) {
516 0 : if (!pSpacePoints->Check(clid->fId)) {
517 0 : HLTError("can not find spacepoint 0x%08x", clid->fId);
518 : continue;
519 : }
520 0 : if (writtenClusterIds) {
521 0 : writtenClusterIds->push_back(clid->fId);
522 0 : }
523 :
524 : // FIXME: there is a bug in the calculation of the residuals stored with the
525 : // assiciated space point, calculate again, but needs to be fixed
526 0 : float deltapad =pSpacePoints->GetY(clid->fId)-clrow->GetU();//clid->fdU;
527 0 : float deltatime =pSpacePoints->GetZ(clid->fId)-clrow->GetV();//clid->fdV;
528 0 : if (TMath::Abs(deltapad)>=AliHLTTPCDefinitions::fgkMaxClusterDeltaPad ||
529 0 : TMath::Abs(deltatime)>=AliHLTTPCDefinitions::fgkMaxClusterDeltaTime) {
530 0 : AliHLTUInt8_t slice = AliHLTTPCSpacePointData::GetSlice(clid->fId);
531 0 : AliHLTUInt8_t partition = AliHLTTPCSpacePointData::GetPatch(clid->fId);
532 0 : HLTFatal("cluster 0x%08x slice %d partition %d: residual out of range - pad %f (max %d), time %f (max %d)", clid->fId, slice, partition, deltapad, AliHLTTPCDefinitions::fgkMaxClusterDeltaPad, deltatime, AliHLTTPCDefinitions::fgkMaxClusterDeltaTime);
533 0 : }
534 0 : float sigmaY2=pSpacePoints->GetYWidth(clid->fId);
535 0 : float sigmaZ2=pSpacePoints->GetZWidth(clid->fId);
536 0 : AliHLTUInt64_t charge=(AliHLTUInt64_t)pSpacePoints->GetCharge(clid->fId);
537 0 : AliHLTUInt64_t qmax=(AliHLTUInt64_t)pTPCRawSpacePoints->GetQMax(clid->fId);
538 : // cout << " row " << setfill(' ') << setw(3) << fixed << right << row
539 : // << " pad " << setfill(' ') << setw(7) << fixed << right << setprecision (4) << pSpacePoints->GetY(clid->fId)
540 : // << " dpad " << setfill(' ') << setw(7) << fixed << right << setprecision (4) << deltapad
541 : // << " time " << setfill(' ') << setw(7) << fixed << right << setprecision (4) << pSpacePoints->GetZ(clid->fId)
542 : // << " dtime " << setfill(' ') << setw(7) << fixed << right << setprecision (4) << deltatime
543 : // << " charge " << setfill(' ') << setw(5) << fixed << right << charge
544 : // << " qmax " << setfill(' ') << setw(4) << fixed << right << qmax
545 : // << endl;
546 :
547 : // time and pad coordinates are scaled and transformed to integer values for
548 : // both cluster and track point before calculating the residual. this makes
549 : // the compression lossless with respect to the format without track model
550 : // compression
551 0 : AliHLTUInt64_t deltapad64=0;
552 0 : AliHLTUInt32_t signDeltaPad=0;
553 0 : if (!isnan(deltapad)) {
554 0 : double clusterpad=pSpacePoints->GetY(clid->fId);
555 0 : double trackpad=clrow->GetU();
556 0 : if (clusterpad<0.) {
557 0 : HLTError("cluster 0x%08x has negative pad position", clid->fId, clusterpad);
558 : }
559 0 : clusterpad*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualPad].fScale;
560 0 : trackpad*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualPad].fScale;
561 0 : AliHLTUInt64_t clusterpad64=(AliHLTUInt64_t)round(clusterpad);
562 : AliHLTUInt64_t trackpad64=0;
563 0 : if (trackpad>0.) trackpad64=(AliHLTUInt64_t)round(trackpad);
564 0 : if (clusterpad64<trackpad64) {
565 0 : deltapad64=trackpad64-clusterpad64;
566 0 : signDeltaPad=1;
567 0 : } else {
568 0 : deltapad64=clusterpad64-trackpad64;
569 0 : signDeltaPad=0;
570 : }
571 0 : }
572 0 : AliHLTUInt64_t deltatime64=0;
573 0 : AliHLTUInt32_t signDeltaTime=0;
574 0 : if (!isnan(deltatime)) {
575 0 : double clustertime=pSpacePoints->GetZ(clid->fId);
576 0 : double tracktime=clrow->GetV();
577 0 : if (clustertime<0.) {
578 0 : HLTError("cluster 0x%08x has negative time position", clid->fId, clustertime);
579 : }
580 0 : clustertime*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualTime].fScale;
581 0 : tracktime*=AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kResidualTime].fScale;
582 0 : AliHLTUInt64_t clustertime64=(AliHLTUInt64_t)round(clustertime);
583 : AliHLTUInt64_t tracktime64=0;
584 0 : if (tracktime>0.) tracktime64=(AliHLTUInt64_t)round(tracktime);
585 0 : if (clustertime64<tracktime64) {
586 0 : deltatime64=tracktime64-clustertime64;
587 0 : signDeltaTime=1;
588 0 : } else {
589 0 : deltatime64=clustertime64-tracktime64;
590 0 : signDeltaTime=0;
591 : }
592 0 : }
593 0 : AliHLTUInt64_t sigmaY264=0;
594 0 : if (!isnan(sigmaY2)) sigmaY264=(AliHLTUInt64_t)round(sigmaY2*AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaY2].fScale);
595 : // we can safely use the upper limit as this is an unphysical cluster, no impact to physics
596 0 : if (sigmaY264 >= (unsigned)1<<AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaY2].fBitLength) {
597 0 : sigmaY264 = (1<<AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaY2].fBitLength)-1;
598 0 : }
599 0 : AliHLTUInt64_t sigmaZ264=0;
600 0 : if (!isnan(sigmaZ2)) sigmaZ264=(AliHLTUInt64_t)round(sigmaZ2*AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaZ2].fScale);
601 : // we can safely use the upper limit as this is an unphysical cluster, no impact to physics
602 0 : if (sigmaZ264 >= (unsigned)1<<AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaZ2].fBitLength) {
603 0 : sigmaZ264 = (1<<AliHLTTPCDefinitions::fgkClusterParameterDefinitions[AliHLTTPCDefinitions::kSigmaZ2].fBitLength)-1;
604 0 : }
605 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualPad , deltapad64);
606 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaPad);
607 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kResidualTime , deltatime64);
608 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(signDeltaTime);
609 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kSigmaY2 , sigmaY264);
610 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kSigmaZ2 , sigmaZ264);
611 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kCharge , charge);
612 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputParameterBits(AliHLTTPCDefinitions::kQMax , qmax);
613 0 : if (bWriteSuccess) writtenClusters++;
614 0 : }
615 0 : }
616 :
617 : // set to next trackpoint
618 0 : if (bReverse) {
619 0 : if (clrow!=fRawTrackPoints.begin()) {
620 0 : AliHLTUInt32_t nextPartition=AliHLTTPCSpacePointData::GetPatch((clrow-1)->GetId());
621 0 : AliHLTUInt32_t nextTrackRow=AliHLTTPCSpacePointData::GetNumber((clrow-1)->GetId());
622 0 : nextTrackRow+=AliHLTTPCGeometry::GetFirstRow(nextPartition);
623 0 : if (thisTrackRow+1==nextTrackRow) {
624 0 : clrow--;
625 0 : } else {
626 : // switch direction start from beginning
627 0 : clrow=fRawTrackPoints.begin();
628 : bReverse=false;
629 : }
630 0 : } else {
631 : // all trackpoints processed
632 0 : clrow=fRawTrackPoints.end();
633 : }
634 : } else {
635 0 : clrow++;
636 : }
637 : continue;
638 0 : } else {
639 : // sequence not ordered, search
640 : // this has been fixed and the search is no longer necessary
641 : // for (clrow=fRawTrackPoints.begin(); clrow!=fRawTrackPoints.end(); clrow++) {
642 : // if ((AliHLTTPCSpacePointData::GetNumber(clrow->GetId())+AliHLTTPCGeometry::GetFirstRow(AliHLTTPCSpacePointData::GetPatch(clrow->GetId())))==row) break;
643 : // }
644 : // if (clrow==fRawTrackPoints.end()) {
645 : // clrow=fRawTrackPoints.begin();
646 : // HLTWarning("no trackpoint on row %d, current point %d", row, thisTrackRow);
647 : // }
648 : }
649 0 : }
650 : // no cluster on that padrow
651 0 : AliHLTUInt32_t haveClusters=0;
652 0 : bWriteSuccess=bWriteSuccess && pDeflater->OutputBit(haveClusters);
653 0 : }
654 :
655 0 : if (!bWriteSuccess) {
656 : // TODO: code review 2015-02-10 Matthias.Richter@scieq.net
657 : // misleading error code, there are two reasons for failed write operation
658 : // 1) target buffer overflow -> -ENOSPC
659 : // 2) value range excess -> -ERANGE
660 0 : return -ENOSPC;
661 : }
662 :
663 : int allClusters=0;
664 0 : for (clrow=fRawTrackPoints.begin(); clrow!=fRawTrackPoints.end(); clrow++) {
665 0 : allClusters+=clrow->GetSpacepoints().size();
666 : }
667 0 : if (allClusters!=writtenClusters) {
668 0 : HLTError("track %d mismatch in written clusters: %d but expected %d", GetTrackId(), writtenClusters, allClusters);
669 : }
670 :
671 0 : pDeflater->Pad8Bits();
672 0 : return pDeflater->GetCurrentByteOutputPosition()-dataPosition;
673 0 : }
674 :
675 : int AliHLTTPCTrackGeometry::Read(const AliHLTUInt8_t* buffer,
676 : AliHLTUInt32_t size,
677 : float bz,
678 : AliHLTUInt32_t& clusterBlockSize,
679 : const char* /*option*/)
680 : {
681 : // read track block from buffer
682 : int iResult=0;
683 0 : if (!buffer) return -EINVAL;
684 0 : if (size<sizeof(AliHLTTPCTrackBlock)) {
685 0 : HLTError("buffer does not contain valid data of track model clusters");
686 0 : return -ENODATA;
687 : }
688 0 : const AliHLTTPCTrackBlock* pTrackBlock=reinterpret_cast<const AliHLTTPCTrackBlock*>(buffer);
689 0 : if (pTrackBlock->fSize>size) {
690 0 : HLTError("inconsistent track data block of size %d exceeds available buffer of size %d", pTrackBlock->fSize, size);
691 0 : return -ENODATA;
692 : }
693 0 : if (pTrackBlock->fSize<sizeof(AliHLTTPCTrackBlock)) {
694 0 : HLTError("inconsistent size of track data block specified in the header: %d", pTrackBlock->fSize);
695 0 : return -ENODATA;
696 : }
697 0 : AliHLTExternalTrackParam param;
698 0 : memset(¶m, 0, sizeof(param));
699 0 : param.fAlpha =( pTrackBlock->fSlice + 0.5 ) * TMath::Pi() / 9.0;
700 0 : if (param.fAlpha>TMath::TwoPi()) param.fAlpha-=TMath::TwoPi();
701 0 : param.fX = pTrackBlock->fX;
702 0 : param.fY = pTrackBlock->fY;
703 0 : param.fZ = pTrackBlock->fZ;
704 0 : param.fSinPsi = pTrackBlock->fSinPsi;
705 0 : param.fTgl = pTrackBlock->fTgl;
706 0 : param.fq1Pt = pTrackBlock->fq1Pt;
707 0 : AliHLTGlobalBarrelTrack track(param);
708 0 : if ((iResult=track.CalculateHelixParams(bz))<0) {
709 0 : HLTError("failed to calculate helix params: %d", iResult);
710 0 : return iResult;
711 : }
712 0 : if ((iResult=CalculateTrackPoints(track))<0) {
713 0 : HLTError("failed to calculate track points: %d", iResult);
714 0 : return iResult;
715 : }
716 0 : clusterBlockSize=pTrackBlock->fSize-sizeof(AliHLTTPCTrackBlock);
717 0 : return sizeof(AliHLTTPCTrackBlock);
718 0 : }
|
