3 //**************************************************************************
4 //* This file is property of and copyright by the ALICE HLT Project *
5 //* ALICE Experiment at CERN, All rights reserved. *
7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
8 //* for The ALICE HLT Project. *
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 //**************************************************************************
19 /** @file AliHLTGlobalBarrelTrack.cxx
20 @author Matthias Richter
22 @brief An AliKalmanTrack implementation for global HLT barrel tracks.
25 // see header file for class documentation
27 // refer to README to build package
29 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
34 #include "AliHLTGlobalBarrelTrack.h"
35 #include "AliHLTSpacePointContainer.h"
36 #include "AliHLTTrackGeometry.h"
37 #include "AliHLTMisc.h"
38 #include "TClonesArray.h"
45 /** ROOT macro for the implementation of ROOT specific class methods */
46 ClassImp(AliHLTGlobalBarrelTrack)
48 AliHLTGlobalBarrelTrack::AliHLTGlobalBarrelTrack()
60 // see header file for class documentation
62 // refer to README to build package
64 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
67 AliHLTGlobalBarrelTrack::AliHLTGlobalBarrelTrack(const AliHLTGlobalBarrelTrack& t)
70 , fLastX(t.GetLastPointX())
71 , fLastY(t.GetLastPointY())
72 , fTrackID(t.TrackID())
73 , fHelixRadius(t.fHelixRadius)
74 , fHelixCenterX(t.fHelixCenterX)
75 , fHelixCenterY(t.fHelixCenterY)
79 // see header file for class documentation
80 fPoints.assign(t.fPoints.begin(), t.fPoints.end());
83 AliHLTGlobalBarrelTrack::AliHLTGlobalBarrelTrack(const AliHLTExternalTrackParam& p )
88 , fTrackID(p.fTrackID)
95 // see header file for class documentation
97 // the 5 track parameters are named in the AliHLTExternalTrackParam
98 // while AliExternalTrackParam just uses an array[5]
99 // the members have the same order, fY is the first one
100 Set(p.fX, p.fAlpha, &p.fY, p.fC);
101 SetPoints(p.fPointIDs, p.fNPoints);
102 SetNumberOfClusters(p.fNPoints);
103 //SetIntegratedLength(GetPathLengthTo( GetLastPointX(), b);
106 AliHLTGlobalBarrelTrack::AliHLTGlobalBarrelTrack(const AliExternalTrackParam& p )
118 // see header file for class documentation
119 *(dynamic_cast<AliExternalTrackParam*>(this))=p;
122 AliHLTGlobalBarrelTrack::~AliHLTGlobalBarrelTrack()
124 // see header file for class documentation
128 Double_t AliHLTGlobalBarrelTrack::GetPathLengthTo( Double_t x, Double_t b ) const
130 // calculate the trajectory length for dx step
132 Double_t dx = x - GetX();
133 Double_t ey = GetSnp();
134 if( TMath::Abs( ey )>=kAlmost1 ) return 0;
136 Double_t ex = TMath::Sqrt(1-ey*ey);
137 Double_t k = GetC(b);
139 Double_t ey1 = k * dx + ey;
141 // check for intersection with X=x
143 if ( TMath::Abs( ey1 ) >= kAlmost1 ) return 0;
145 Double_t ex1 = TMath::Sqrt(1-ey1*ey1);
147 Double_t ss = ey + ey1;
148 Double_t cc = ex + ex1;
150 if ( TMath::Abs( cc ) < 1.e-4 ) return 0;
152 Double_t tg = ss / cc;
153 Double_t dl = dx * TMath::Sqrt( 1 + tg * tg );
154 Double_t dSin = dl * k / 2;
155 if ( dSin > 1 ) dSin = 1;
156 if ( dSin < -1 ) dSin = -1;
157 Double_t dS = ( TMath::Abs( k ) > 1.e-4 ) ? ( 2 * TMath::ASin( dSin ) / k ) : dl;
159 return dS*TMath::Sqrt(1 + GetTgl()*GetTgl() );
165 int AliHLTGlobalBarrelTrack::ConvertTrackDataArray(const AliHLTTracksData* pTracks, unsigned sizeInByte, vector<AliHLTGlobalBarrelTrack> &tgtArray )
167 // see header file for class documentation
169 if (!pTracks || sizeInByte<sizeof(AliHLTTracksData) || pTracks->fCount==0) return 0;
171 const AliHLTUInt8_t* pEnd=reinterpret_cast<const AliHLTUInt8_t*>(pTracks);
174 tgtArray.resize(pTracks->fCount + tgtArray.size());
175 const AliHLTUInt8_t* pCurrent=reinterpret_cast<const AliHLTUInt8_t*>(pTracks->fTracklets);
176 for (unsigned i=0; i<pTracks->fCount; i++) {
177 if (pCurrent+sizeof(AliHLTExternalTrackParam)>pEnd) {
178 iResult=-EINVAL; break;
180 const AliHLTExternalTrackParam* track=reinterpret_cast<const AliHLTExternalTrackParam*>(pCurrent);
181 if (pCurrent+sizeof(AliHLTExternalTrackParam)+track->fNPoints*sizeof(UInt_t)>pEnd) {
182 iResult=-EINVAL; break;
185 pCurrent+=sizeof(AliHLTExternalTrackParam)+track->fNPoints*sizeof(UInt_t);
187 if (iResult<0) tgtArray.clear();
188 else iResult=tgtArray.size();
192 UInt_t AliHLTGlobalBarrelTrack::GetNumberOfPoints() const
194 // see header file for class documentation
195 return fPoints.size();
198 const UInt_t* AliHLTGlobalBarrelTrack::GetPoints() const
200 // see header file for class documentation
201 if (fPoints.size()==0) return NULL;
205 int AliHLTGlobalBarrelTrack::SetPoints(const UInt_t* pArray, UInt_t arraySize)
207 // see header file for class documentation
208 if (!pArray || arraySize==0) return 0;
209 fPoints.resize(arraySize);
210 for (unsigned i=0; i<arraySize; i++) fPoints[i]=pArray[i];
211 return fPoints.size();
214 int AliHLTGlobalBarrelTrack::CalculateHelixParams()
216 // calculate radius and center of the helix
217 float bfield=AliHLTMisc::Instance().GetBz();
218 if (TMath::Abs(bfield)<kAlmost0) {
219 // no magnetic field -> straight lines
220 fHelixRadius=kVeryBig;
221 fHelixCenterX=kVeryBig;
222 fHelixCenterY=kVeryBig;
224 fHelixRadius = GetSignedPt()/(-kB2C*bfield);
225 Double_t trackPhi = Phi()-GetAlpha();
227 //cout << "Helix: phi=" << trackPhi << " x=" << GetX() << " y=" << GetY() << endl;
228 fHelixCenterX = GetX() + fHelixRadius * sin(trackPhi);
229 fHelixCenterY = GetY() - fHelixRadius * cos(trackPhi);
230 //cout << "Helix: center" << " x=" << fHelixCenterX << " y=" << fHelixCenterY << endl;
235 int AliHLTGlobalBarrelTrack::CalculateCrossingPoint(float xPlane, float alphaPlane, float& y, float& z)
237 // calculate crossing point of helix with a plane in yz
238 // in the local coordinates of the plane
240 if (TMath::Abs(fHelixRadius)<kAlmost0 &&
241 (iResult=CalculateHelixParams())<0) {
245 if (TMath::Abs(fHelixRadius)>=kVeryBig) {
246 // no magnetic field -> straight lines
248 // rotate helix center to local coordinates of the plane reference frame
249 float cosa=TMath::Cos(alphaPlane-GetAlpha());
250 float sina=TMath::Sin(alphaPlane-GetAlpha());
251 float cx= fHelixCenterX * cosa + fHelixCenterY * sina;
252 float cy=-fHelixCenterX * sina + fHelixCenterY * cosa;
254 // crossing point of helix with plane
255 Double_t aa = (cx - xPlane)*(cx - xPlane);
256 Double_t r2 = fHelixRadius*fHelixRadius;
257 if(aa > r2) // no crossing
260 Double_t aa2 = sqrt(r2 - aa);
261 Double_t y1 = cy + aa2;
262 Double_t y2 = cy - aa2;
264 if(TMath::Abs(y2) < TMath::Abs(y1)) y = y2;
266 // calculate the arc length between (x,y) and (x0,y0) with radius (cx,cy)
267 // reference point is (x0,y0) rotated by the diffence of the plane angle and
268 // the track reference frame angle alpha
270 Double_t angle1 = atan2((y - cy),(xPlane - cx));
271 if(angle1 < 0) angle1 += TMath::TwoPi();
273 // 2) angle of (x0,y0)
274 float x0= GetX() * cosa + GetY() * sina;
275 float y0=-GetX() * sina + GetY() * cosa;
276 Double_t angle2 = atan2((y0 - cy),(x0 - cx));
277 if(angle2 < 0) angle2 += TMath::TwoPi();
279 // 3) angle between (x,y) and (x0,y0)
280 Double_t diffangle = angle1 - angle2;
281 diffangle = fmod(diffangle,TMath::TwoPi());
284 Double_t arclength = TMath::Abs(diffangle*fHelixRadius);
286 // 5) direction depending on whether going outwards or inwards
287 int direction=GetX()>xPlane?-1:1;
288 z = GetZ() + direction*arclength*GetTgl();
290 //cout << "x=" << xPlane << " y=" << y << " cx=" << cx << " cy=" << cy << " a1=" << angle1 << " a2=" << angle2 << " diffa=" << diffangle << " s=" << arclength << " z=" << z << endl;
295 void AliHLTGlobalBarrelTrack::Print(Option_t* option) const
297 // see header file for class documentation
298 cout << "********* Track Id: " << fTrackID << " *******************" << endl;
299 AliExternalTrackParam::Print(option);
300 // cout << " Alpha " << GetAlpha();
301 // cout << " X " << GetX();
302 // cout << " Y " << GetY();
303 // cout << " Z " << GetZ() << endl;
304 // cout << " Snp " << GetSnp();
305 // cout << " Tgl " << GetTgl();
306 // cout << " Signed1Pt " << GetSigned1Pt() << endl;
309 void AliHLTGlobalBarrelTrack::Draw(Option_t *option)
311 /// Inherited from TObject, draw the track
313 float center[2]={0.5,0.5};
315 if (TMath::Abs(fHelixRadius)<kAlmost0 &&
316 (CalculateHelixParams())<0) {
320 TString strOption(option);
321 if (strOption.IsNull()) strOption="spacepoints trackarc";
322 std::auto_ptr<TObjArray> tokens(strOption.Tokenize(" "));
323 if (!tokens.get()) return;
324 for (int i=0; i<tokens->GetEntriesFast(); i++) {
325 if (!tokens->At(i)) continue;
327 TString arg=tokens->At(i)->GetName();
330 if (arg.BeginsWith(key)) {
331 arg.ReplaceAll(key, "");
336 if (arg.BeginsWith(key)) {
337 arg.ReplaceAll(key, "");
338 center[0]=arg.Atof();
342 if (arg.BeginsWith(key)) {
343 arg.ReplaceAll(key, "");
344 center[1]=arg.Atof();
348 if (arg.CompareTo(key)==0) {
349 if (fSpacePoints) DrawProjXYSpacePoints(option, fSpacePoints, scale, center);
353 if (arg.CompareTo(key)==0) {
354 DrawProjXYTrack(option, scale, center);
360 int AliHLTGlobalBarrelTrack::DrawProjXYSpacePoints(Option_t */*option*/, const AliHLTSpacePointContainer* spacePoints, const float scale, float center[2])
362 /// draw space points
365 if (!spacePoints) return -EINVAL;
367 const UInt_t* pointids=GetPoints();
368 for (unsigned i=0; i<GetNumberOfPoints() && pointids; i++) {
369 float clusterphi=spacePoints->GetPhi(pointids[i]);
370 float cosphi=TMath::Cos(clusterphi);
371 float sinphi=TMath::Sin(clusterphi);
372 float clusterx=spacePoints->GetX(pointids[i]);
373 float clustery=spacePoints->GetY(pointids[i]);
375 float pointx= clusterx*sinphi + clustery*cosphi;
376 float pointy=-clusterx*cosphi + clustery*sinphi;
378 // FIXME: cleanup of marker objects
379 TMarker* m=new TMarker(pointx/(2*scale)+center[0], pointy/(2*scale)+center[1], 3);
380 m->SetMarkerColor(markerColor);
386 // FIXME: make this a general geometry definition
387 // through an abstract class interface
388 const Double_t gkTPCX[159] = {
550 int AliHLTGlobalBarrelTrack::DrawProjXYTrack(Option_t *option, const float scale, float center[2])
553 bool bDrawArc=false; // draw TArc
554 bool bNoTrackPoints=false; // don't draw track points
555 TString strOption(option);
556 if (strOption.IsNull()) strOption="spacepoints trackarc";
557 std::auto_ptr<TObjArray> tokens(strOption.Tokenize(" "));
558 if (!tokens.get()) return 0;
559 for (int i=0; i<tokens->GetEntriesFast(); i++) {
560 if (!tokens->At(i)) continue;
562 TString arg=tokens->At(i)->GetName();
565 if (arg.BeginsWith(key)) {
570 if (arg.BeginsWith(key)) {
576 float cosa=TMath::Cos(GetAlpha());
577 float sina=TMath::Sin(GetAlpha());
581 firstpoint[0]= GetX()*sina + GetY()*cosa;
582 firstpoint[1]=-GetX()*cosa + GetY()*sina;
584 //cout << " first point alpha=" << GetAlpha() << " x: " << firstpoint[0] << " y: " << firstpoint[1] << endl;
585 // FIXME: cleanup of marker objects
586 TMarker* m=new TMarker(firstpoint[0]/(2*scale)+center[0], firstpoint[1]/(2*scale)+center[1], 29);
588 m->SetMarkerColor(2);
592 // draw points in step width and remember the last point
593 float lastpoint[2]={0.0, 0.0};
595 for (; firstpadrow<159 && gkTPCX[firstpadrow]<GetX(); firstpadrow++);
596 DrawProjXYTrackPoints(option, scale, center, firstpadrow, -1, firstpoint);
597 DrawProjXYTrackPoints(option, scale, center, firstpadrow, 1, lastpoint);
600 if (TMath::Abs(fHelixRadius)>=kVeryBig) {
601 // no magnetic field -> straight lines
603 // rotate helix center to local coordinates of the plane reference frame
604 float cx= fHelixCenterX * sina + fHelixCenterY * cosa;
605 float cy=-fHelixCenterX * cosa + fHelixCenterY * sina;
607 float diffx=cx-firstpoint[0];
608 float diffy=cy-firstpoint[1];
610 float phimax=2*TMath::Pi();
611 if (TMath::Abs(diffx)<kAlmost0) {
612 phimin=TMath::Pi()/2;
614 phimin=TMath::ATan(diffy/diffx);
616 if (diffx>0) phimin+=TMath::Pi();
618 diffx=cx-lastpoint[0];
619 diffy=cy-lastpoint[1];
620 if (TMath::Abs(diffx)<kAlmost0) {
621 phimax=TMath::Pi()/2;
623 phimax=TMath::ATan(diffy/diffx);
625 //cout << "diffx=" << diffx << " diffy=" << diffy << " phimin=" << phimin << " phimax=" << phimax << endl;
626 if (diffx>0) phimax+=TMath::Pi();
627 if (phimax<0 && phimin>=0 &&
628 phimax+TMath::TwoPi()-phimin<TMath::Pi()) phimax+=TMath::TwoPi();
629 //if (phimax<0 && TMath::Abs(phimax-phimin)<TMath::Pi()) phimax+=TMath::TwoPi();
630 if (phimax-phimin>TMath::Pi()) phimax-=TMath::TwoPi();
632 if (0/*phimin>phimax*/) {
637 phimin*=360.0/(2*TMath::Pi());
638 phimax*=360.0/(2*TMath::Pi());
639 //cout << " cx=" << cx << " cy=" << cy << " r=" << fHelixRadius << " phimin=" << phimin << " phimax=" << phimax << endl;
640 // FIXME: cleanup of graphics objects
641 TArc* tarc=new TArc(cx/(2*scale)+center[0], cy/(2*scale)+center[1], TMath::Abs(fHelixRadius)/(2*scale), phimin, phimax);
643 tarc->SetFillStyle(0);
650 int AliHLTGlobalBarrelTrack::DrawProjXYTrackPoints(Option_t */*option*/, const float scale, const float center[2], int firstpadrow, int step, float point[2])
652 // draw points in step width and return the last point
653 float offsetAlpha=0.0;
654 float cosa=TMath::Cos(GetAlpha());
655 float sina=TMath::Sin(GetAlpha());
657 for (int padrow=firstpadrow; padrow>=0 && padrow<159; padrow+=step) {
658 float x=gkTPCX[padrow];
666 if ((result=CalculateCrossingPoint(x, GetAlpha()-offsetAlpha, y, z))<1) break;
667 float pointAlpha=TMath::ATan(y/x);
668 if (TMath::Abs(pointAlpha)>TMath::Pi()/18) {
669 offsetAlpha+=(pointAlpha>0?-1:1)*TMath::Pi()/9;
672 cosa=TMath::Cos(GetAlpha()-offsetAlpha);
673 sina=TMath::Sin(GetAlpha()-offsetAlpha);
675 } while (result==0 && shift++<maxshift);
676 if (result<1) continue;
677 point[0]= x*sina + y*cosa;
678 point[1]=-x*cosa + y*sina;
680 //cout << x << " : x=" << x << " y=" << y << endl;
681 // FIXME: cleanup of TMarker objects?
682 TMarker* m=new TMarker(point[0]/(2*scale)+center[0], point[1]/(2*scale)+center[1], z>=0?2:5);
683 m->SetMarkerColor(markerColor);
689 void AliHLTGlobalBarrelTrack::SetTrackGeometry(AliHLTTrackGeometry* points)
691 /// set the instance to the track points container
693 if (!fTrackPoints) return;
694 fTrackPoints->SetTrackId(GetID());
697 int AliHLTGlobalBarrelTrack::AssociateSpacePoints(AliHLTTrackGeometry* trackpoints, AliHLTSpacePointContainer& spacepoints) const
699 /// associate the track space points to the calculated track points
700 AliHLTTrackGeometry* instance=trackpoints;
701 if (!instance) instance=fTrackPoints;
702 if (!instance) return 0;
704 UInt_t nofIds=GetNumberOfPoints();
705 const UInt_t* ids=GetPoints();
706 int result=instance->AssociateSpacePoints(ids, nofIds, spacepoints);
710 int AliHLTGlobalBarrelTrack::ReadTracks(const char* filename, TClonesArray& tgt, AliHLTComponentDataType /*dt*/, unsigned /*specification*/)
712 // open block from file and add to collection
713 if (!filename) return -EINVAL;
715 TString input=filename;
716 input+="?filetype=raw";
717 std::auto_ptr<TFile> pFile(new TFile(input));
718 if (!pFile.get()) return -ENOMEM;
719 if (pFile->IsZombie()) return -ENOENT;
723 std::auto_ptr<TArrayC> buffer(new TArrayC);
724 if (!buffer.get()) return -ENOMEM;
726 buffer->Set(pFile->GetSize());
727 if (pFile->ReadBuffer(buffer->GetArray(), buffer->GetSize())==0) {
728 const AliHLTTracksData* pTracks=reinterpret_cast<const AliHLTTracksData*>(buffer->GetArray());
729 vector<AliHLTGlobalBarrelTrack> tracks;
730 iResult=ConvertTrackDataArray(pTracks, buffer->GetSize(), tracks);
732 int offset=tgt.GetEntriesFast();
733 tgt.ExpandCreate(offset+tracks.size());
734 for (unsigned i=0; i<tracks.size(); i++) {
735 new (tgt[offset+i]) AliHLTGlobalBarrelTrack(tracks[i]);
737 iResult=tracks.size();
739 //HLTError("failed to convert tracks from file %s size %d byte(s) ", filename, pFile->GetSize());
742 //HLTError("failed reading %d byte(s) from file %s", pFile->GetSize(), filename);
749 int AliHLTGlobalBarrelTrack::ReadTrackList(const char* listfile, TClonesArray& tgt, AliHLTComponentDataType dt, unsigned specification)
751 // open blank separated list of files and read tracks
752 ifstream list(listfile);
753 if (!list.good()) return -ENOENT;
757 while (file.ReadLine(list)) {
758 //HLTInfo("adding tracks from file %s", file.Data());
759 int iResult=ReadTracks(file.Data(), tgt, dt, specification);
761 //HLTInfo("failed to add data from file %s: error %d", file.Data(), iResult);