c8c81c5f00c2dc774ba7d7ecebf2c40adfefd534
[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking-ca / AliHLTTPCCATrackerComponent.cxx
1 // @(#) $Id$
2 /**************************************************************************
3  * This file is property of and copyright by the ALICE HLT Project        * 
4  * ALICE Experiment at CERN, All rights reserved.                         *
5  *                                                                        *
6  * Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de>        *
7  *                  Ivan Kisel <kisel@kip.uni-heidelberg.de>              *
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 ///////////////////////////////////////////////////////////////////////////////
20 //                                                                           //
21 // a TPC tracker processing component for the HLT based on CA by Ivan Kisel  //
22 //                                                                           //
23 ///////////////////////////////////////////////////////////////////////////////
24
25 #if __GNUC__>= 3
26 using namespace std;
27 #endif
28
29 #include "AliHLTTPCCATrackerComponent.h"
30 #include "AliHLTTPCTransform.h"
31 #include "AliHLTTPCCATracker.h"
32 #include "AliHLTTPCCAHit.h"
33 #include "AliHLTTPCCAOutTrack.h"
34
35 #include "AliHLTTPCSpacePointData.h"
36 #include "AliHLTTPCClusterDataFormat.h"
37 #include "AliHLTTPCTransform.h"
38 #include "AliHLTTPCTrackSegmentData.h"
39 #include "AliHLTTPCTrackArray.h"
40 #include "AliHLTTPCTrackletDataFormat.h"
41 #include "AliHLTTPCDefinitions.h"
42 #include "TMath.h"
43 #include "AliTPC.h"
44 #include "AliTPCParam.h"
45 #include "AliRun.h"
46 #include <stdlib.h>
47 #include <iostream>
48 #include <errno.h>
49
50
51 // this is a global object used for automatic component registration, do not use this
52 AliHLTTPCCATrackerComponent gAliHLTTPCCATrackerComponent;
53
54 ClassImp(AliHLTTPCCATrackerComponent)
55
56 AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent()
57   :
58   fTracker(NULL),
59   fBField(0)
60 {
61   // see header file for class documentation
62   // or
63   // refer to README to build package
64   // or
65   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
66 }
67
68 AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent(const AliHLTTPCCATrackerComponent&)
69   :
70   fTracker(NULL),
71   fBField(0)
72 {
73   // see header file for class documentation
74   HLTFatal("copy constructor untested");
75 }
76
77 AliHLTTPCCATrackerComponent& AliHLTTPCCATrackerComponent::operator=(const AliHLTTPCCATrackerComponent&)
78 {
79   // see header file for class documentation
80   HLTFatal("assignment operator untested");
81   return *this;
82 }
83
84 AliHLTTPCCATrackerComponent::~AliHLTTPCCATrackerComponent()
85 {
86   // see header file for class documentation
87   delete fTracker;
88 }
89
90 //
91 // Public functions to implement AliHLTComponent's interface.
92 // These functions are required for the registration process
93 //
94
95 const char* AliHLTTPCCATrackerComponent::GetComponentID() 
96 {
97   // see header file for class documentation
98   return "TPCCATracker";
99 }
100
101 void AliHLTTPCCATrackerComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list) 
102 {
103   // see header file for class documentation
104   list.clear();
105   list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
106 }
107
108 AliHLTComponentDataType AliHLTTPCCATrackerComponent::GetOutputDataType() 
109 {
110   // see header file for class documentation
111   return AliHLTTPCDefinitions::fgkTrackSegmentsDataType;
112 }
113
114 void AliHLTTPCCATrackerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) 
115 {
116   // define guess for the output data size
117   constBase = 200;       // minimum size
118   inputMultiplier = 0.5; // size relative to input 
119 }
120
121 AliHLTComponent* AliHLTTPCCATrackerComponent::Spawn() 
122 {
123   // see header file for class documentation
124   return new AliHLTTPCCATrackerComponent;
125 }
126
127 int AliHLTTPCCATrackerComponent::DoInit( int argc, const char** argv )
128 {
129   // Initialize the CA tracker component 
130   //
131   // arguments could be:
132   // bfield - the magnetic field value
133   //
134
135   if ( fTracker ) return EINPROGRESS;
136   
137   fTracker = new AliHLTTPCCATracker();
138   
139   // read command line
140
141   int i = 0;
142   char* cpErr;
143   while ( i < argc ){
144     if ( !strcmp( argv[i], "bfield" ) ){
145       if ( i+1 >= argc )
146         {
147           Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing B-field", "Missing B-field specifier." );
148           return ENOTSUP;
149         }
150       fBField = strtod( argv[i+1], &cpErr );
151       if ( *cpErr )
152         {
153           Logging( kHLTLogError, "HLT::TPCCATracker::DoInit", "Missing multiplicity", "Cannot convert B-field specifier '%s'.", argv[i+1] );
154           return EINVAL;
155         }
156
157       Logging( kHLTLogDebug, "HLT::TPCCATracker::DoInit", "Reading command line",
158                "Magnetic field value is set to %f kG", fBField );
159
160       i += 2;
161       continue;
162     }
163     
164     Logging(kHLTLogError, "HLT::TPCCATracker::DoInit", "Unknown Option", "Unknown option '%s'", argv[i] );
165     return EINVAL;
166   }
167   
168   return 0;
169 }
170
171 int AliHLTTPCCATrackerComponent::DoDeinit()
172 {
173   // see header file for class documentation
174   if ( fTracker ) delete fTracker;
175   fTracker = NULL;  
176   return 0;
177 }
178
179 int AliHLTTPCCATrackerComponent::DoEvent
180
181  const AliHLTComponentEventData& evtData, 
182  const AliHLTComponentBlockData* blocks, 
183  AliHLTComponentTriggerData& trigData, 
184  AliHLTUInt8_t* outputPtr, 
185  AliHLTUInt32_t& size, 
186  vector<AliHLTComponentBlockData>& outputBlocks )
187 {
188   // Event reconstruction in one TPC slice with CA Tracker
189
190   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "DoEvent", "DoEvent()" );
191   if ( evtData.fBlockCnt<=0 )
192     {
193       Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "no blocks in event" );
194       return 0;
195     }
196   
197   const AliHLTComponentBlockData* iter = NULL;
198   unsigned long ndx;
199   AliHLTTPCClusterData* inPtrSP; 
200  
201   // Determine the slice number 
202   
203   Int_t slice=-1;
204   {
205     std::vector<Int_t> slices;
206     std::vector<Int_t>::iterator slIter;
207     std::vector<unsigned> sliceCnts;
208     std::vector<unsigned>::iterator slCntIter;
209   
210     for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ){
211       iter = blocks+ndx;
212       if ( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
213
214       slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
215       Bool_t found = 0;
216       slCntIter = sliceCnts.begin();
217       for( slIter = slices.begin(); slIter!=slices.end(); slIter++, slCntIter++ ){
218         if ( *slIter == slice ){
219           found = kTRUE;
220           break;
221         }
222       }
223       if ( !found ){
224         slices.push_back( slice );
225         sliceCnts.push_back( 1 );
226       } else *slCntIter++;     
227     }
228   
229     
230     // Determine slice number to really use.
231     if ( slices.size()>1 )
232       {
233         Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
234                  "Multiple slice numbers found in event 0x%08lX (%lu). Determining maximum occuring slice number...",
235                  evtData.fEventID, evtData.fEventID );
236         unsigned maxCntSlice=0;
237         slCntIter = sliceCnts.begin();
238         for( slIter = slices.begin(); slIter != slices.end(); slIter++, slCntIter++ )
239           {
240             Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
241                      "Slice %lu found %lu times.", *slIter, *slCntIter );
242             if ( maxCntSlice<*slCntIter )
243               {
244                 maxCntSlice = *slCntIter;
245                 slice = *slIter;
246               }
247           }
248         Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
249                  "Using slice %lu.", slice );
250       }
251     else if ( slices.size()>0 )
252       {
253         slice = *(slices.begin());
254       }
255   }
256   
257   if( slice<0 ){
258     Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "no slices found in event" );
259     return 0;
260   }
261
262
263   // Initialize the tracker
264
265   Double_t Bz = fBField;
266   
267   {
268     if( !fTracker ) fTracker = new AliHLTTPCCATracker;
269     Int_t iSec = slice;
270     Double_t inRmin = 83.65; 
271     Double_t inRmax = 133.3;
272     Double_t outRmin = 133.5; 
273     Double_t outRmax = 247.7;
274     Double_t plusZmin = 0.0529937; 
275     Double_t plusZmax = 249.778; 
276     Double_t minusZmin = -249.645; 
277     Double_t minusZmax = -0.0799937; 
278     Double_t dalpha = 0.349066;
279     Double_t alpha = 0.174533 + dalpha*iSec;
280     
281     Bool_t zPlus = (iSec<18|| (iSec>=36&&iSec<54) );
282     Bool_t rInner = (iSec<36);
283     Double_t zMin =  zPlus ?plusZmin :minusZmin;
284     Double_t zMax =  zPlus ?plusZmax :minusZmax;
285     Double_t rMin =  rInner ?inRmin :outRmin;
286     Double_t rMax =  rInner ?inRmax :outRmax;
287     Int_t inNRows = 63;
288     Int_t outNRows = 96;
289     Double_t inRowXFirst = 85.225;
290     Double_t outRowXFirst =135.1;
291     Double_t inRowXStep = 0.75;
292     Double_t outRowXStep = 1.;
293     Int_t nRows = rInner ?inNRows :outNRows;
294     Double_t rowXFirst = rInner ?inRowXFirst :outRowXFirst;
295     Double_t rowXStep = rInner ?inRowXStep :outRowXStep;
296     
297     Int_t nSectors = 72/2;
298     
299     Double_t padPitch = 0.4;
300     Double_t sigmaZ = 0.228808;
301     
302     //TPCZmin = -249.645, ZMax = 249.778    
303      
304     AliHLTTPCCAParam param;
305     param.Initialize( iSec, inNRows+outNRows, inRowXFirst, inRowXStep,alpha, dalpha,
306                       inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, Bz );
307       
308     fTracker->Initialize( param );
309
310     // recalculate outer rows -> necessary for a moment
311
312     for( Int_t irow=0; irow<outNRows; irow++){
313       fTracker->Rows()[inNRows+irow].X() = outRowXFirst + irow*outRowXStep;
314     }     
315   }
316
317     
318   // min and max patch numbers and row numbers
319
320   Int_t row[2] = {0,0};
321   Int_t minPatch=INT_MAX, maxPatch = 0;
322
323   // total n Hits
324
325   Int_t nHitsTotal = 0;
326
327   // sort patches
328
329   std::vector<unsigned long> patchIndices;
330
331   for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ){
332     iter = blocks+ndx;      
333     if( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
334     if( slice!=AliHLTTPCDefinitions::GetMinSliceNr( *iter ) ) continue;
335     inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
336     nHitsTotal+=inPtrSP->fSpacePointCnt;
337     Int_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
338     if ( minPatch>patch ){
339       minPatch = patch;
340       row[0] = AliHLTTPCTransform::GetFirstRow( patch );
341     }
342     if ( maxPatch<patch ){
343       maxPatch = patch;
344       row[1] = AliHLTTPCTransform::GetLastRow( patch );
345     }
346     std::vector<unsigned long>::iterator pIter = patchIndices.begin(); 
347     while( pIter!=patchIndices.end() && AliHLTTPCDefinitions::GetMinPatchNr( blocks[*pIter] ) < patch ){
348       pIter++;
349     }
350     patchIndices.insert( pIter, ndx );
351   }
352            
353
354   // pass event to CA Tracker
355   
356   fTracker->StartEvent();
357
358   AliHLTTPCCAHit * vHits = new AliHLTTPCCAHit[nHitsTotal]; // CA hit array
359   Double_t * vHitStoreX = new Double_t[nHitsTotal];       // hit X coordinates
360   Int_t * vHitStoreID = new Int_t[nHitsTotal];            // hit ID's
361
362   Int_t nHits = 0;
363  
364   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
365            "Total %d hits to read for slice %d", nHitsTotal, slice );
366
367   for( std::vector<unsigned long>::iterator pIter = patchIndices.begin(); pIter!=patchIndices.end(); pIter++ ){
368     ndx = *pIter;
369     iter = blocks+ndx;
370       
371     Int_t patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
372     inPtrSP = (AliHLTTPCClusterData*)(iter->fPtr);
373       
374     Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
375              "Reading %d hits for slice %d - patch %d", inPtrSP->fSpacePointCnt, slice, patch );
376       
377     // Read patch hits, row by row
378
379     Int_t oldRow = -1;
380     Int_t nRowHits = 0;
381     Int_t firstRowHit = 0;
382     for (UInt_t i=0; i<inPtrSP->fSpacePointCnt; i++ ){  
383       AliHLTTPCSpacePointData* pSP = &(inPtrSP->fSpacePoints[i]);
384
385       //Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits", "hit pad %d, xyz=%f,%f,%f, sy=%f, sz=%f,", pSP->fPadRow, pSP->fX, pSP->fY, pSP->fZ, TMath::Sqrt(pSP->fSigmaY2), TMath::Sqrt(pSP->fSigmaZ2) );
386
387       if( pSP->fPadRow != oldRow ){
388         if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
389         oldRow = pSP->fPadRow;
390         firstRowHit = nHits;
391         nRowHits = 0;
392       }
393       AliHLTTPCCAHit &h = vHits[nHits];
394       h.Y() = pSP->fY;
395       h.Z() = pSP->fZ;
396       h.ErrY() = TMath::Sqrt(pSP->fSigmaY2);
397       h.ErrZ() = TMath::Sqrt(pSP->fSigmaZ2);  
398       h.ID() = nHits;
399       vHitStoreX[nHits] = pSP->fX;
400       vHitStoreID[nHits] = pSP->fID;
401       nHits++;  
402       nRowHits++;
403     }   
404     if( oldRow>=0 ) fTracker->ReadHitRow( oldRow, vHits+firstRowHit, nRowHits );
405   }
406
407   
408   // reconstruct the event  
409
410   fTracker->Reconstruct();
411
412   Int_t ret = 0;
413
414   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reconstruct",
415            "%d tracks found for slice %d",fTracker->NOutTracks(), slice);
416
417   // write reconstructed tracks
418
419   AliHLTTPCTrackletData* outPtr = (AliHLTTPCTrackletData*)(outputPtr);
420
421   AliHLTTPCTrackSegmentData* currOutTracklet = outPtr->fTracklets;
422
423   Int_t ntracks = fTracker->NOutTracks();
424
425   UInt_t mySize =   ((AliHLTUInt8_t *)currOutTracklet) -  ((AliHLTUInt8_t *)outputPtr);
426
427   outPtr->fTrackletCnt = 0; 
428
429   for( int itr=0; itr<ntracks; itr++ ){
430     
431     AliHLTTPCCAOutTrack &t = fTracker->OutTracks()[itr];    
432
433     //Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Wrtite output","track %d with %d hits", itr, t.NHits());
434
435     // calculate output track size
436
437     UInt_t dSize = sizeof(AliHLTTPCTrackSegmentData) + t.NHits()*sizeof(UInt_t);
438     
439     if( mySize + dSize >size ){
440       Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "Wrtite output","Output buffer size exceed (buffer size %d, current size %d), %d tracks are not stored", size, mySize, ntracks-itr+1);
441       ret = -ENOSPC;
442       break;
443     }
444     
445     // convert CA track parameters to HLT Track Segment
446
447     Int_t iFirstHit = fTracker->OutTrackHits()[t.FirstHitRef()];
448     Int_t iLastHit = fTracker->OutTrackHits()[t.FirstHitRef()+t.NHits()-1];
449     
450     AliHLTTPCCAHit &firstHit = vHits[iFirstHit];
451     AliHLTTPCCAHit &lastHit = vHits[iLastHit];
452    
453     t.Param().TransportBz(Bz, vHitStoreX[iFirstHit], firstHit.Y(), firstHit.Z() );
454
455     currOutTracklet->fX = t.Param().Par()[0];
456     currOutTracklet->fY = t.Param().Par()[1];
457     currOutTracklet->fZ = t.Param().Par()[2];
458     Double_t qp = t.Param().Par()[6];
459     Double_t p = TMath::Abs(qp)>1.e-5 ?1./TMath::Abs(qp) :1.e5;
460     Double_t ex = t.Param().Par()[3];
461     Double_t ey = t.Param().Par()[4];
462     Double_t ez = t.Param().Par()[5];
463     Double_t et = TMath::Sqrt( ex*ex + ey*ey );
464     
465     currOutTracklet->fCharge = (qp>0) ?+1 :(qp<0 ?-1 :0);
466     currOutTracklet->fPt = p*et;
467
468     Double_t h3 =  TMath::Abs(ex) >1.e-5 ? p*ex/et :0;
469     Double_t h4 =  TMath::Abs(ey) >1.e-5 ? p*ey/et :0;
470     Double_t h5;
471     Double_t h6 =  - currOutTracklet->fCharge * p * currOutTracklet->fPt;
472   
473     currOutTracklet->fPterr = ( h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + h6*h6*t.Param().Cov()[27] 
474                                 + 2.*(h3*h4*t.Param().Cov()[13]+h3*h6*t.Param().Cov()[24]+h4*h6*t.Param().Cov()[25] )
475                                 );    
476    currOutTracklet->fPsi = TMath::ATan2(ey, ex);
477     
478     h3 =  ex/(et*et);
479     h4 = -ey/(et*et);
480     currOutTracklet->fPsierr = h3*h3*t.Param().Cov()[9] + h4*h4*t.Param().Cov()[14] + 2.*h3*h4*t.Param().Cov()[13];
481     
482     currOutTracklet->fTgl = TMath::Abs(ex)>1.e-5  ? ez/ex :1.e5;
483       
484     h3 = (TMath::Abs(ex) >1.e-5) ? -ez/ex/ex :0;
485     h5 = (TMath::Abs(ex) >1.e-5) ? 1./ex :0;
486     currOutTracklet->fTglerr =  h3*h3*t.Param().Cov()[9] + h5*h5*t.Param().Cov()[20] + 2.*h3*h5*t.Param().Cov()[18]; 
487     
488     currOutTracklet->fCharge = -currOutTracklet->fCharge;
489     
490     t.Param().TransportBz(Bz, vHitStoreX[iLastHit], lastHit.Y(), lastHit.Z() );
491      
492     currOutTracklet->fLastX = t.Param().Par()[0];
493     currOutTracklet->fLastY = t.Param().Par()[1];
494     currOutTracklet->fLastZ = t.Param().Par()[2];
495
496     currOutTracklet->fNPoints = t.NHits();
497
498     for( Int_t i=0; i<t.NHits(); i++ ){
499       currOutTracklet->fPointIDs[i] = vHitStoreID[fTracker->OutTrackHits()[t.FirstHitRef()+i]];
500     }
501
502     currOutTracklet = (AliHLTTPCTrackSegmentData*)( (Byte_t *)currOutTracklet + dSize );
503     mySize+=dSize;
504     outPtr->fTrackletCnt++; 
505   }
506     
507   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Tracks",
508            "Input: Number of tracks: %d Slice/MinPatch/MaxPatch/RowMin/RowMax: %d/%d/%d/%d/%d.", 
509            ntracks, slice, minPatch, maxPatch, row[0], row[1] );
510   
511   AliHLTComponentBlockData bd;
512   FillBlockData( bd );
513   bd.fOffset = 0;
514   bd.fSize = mySize;
515   bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, minPatch, maxPatch );      
516   outputBlocks.push_back( bd );
517   
518   size = mySize;
519
520   delete [] vHits;
521   delete [] vHitStoreX;
522   delete [] vHitStoreID;
523   
524   return ret;
525
526 }
527
528