1a2e03a23ce2919f0b49b01f3f23fa60ea045f03
[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: Sergey Gorbunov <sergey.gorbunov@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 ///////////////////////////////////////////////////////////////////////////////
22 //                                                                           //
23 // a TPC tracker processing component for the HLT based on CA by Ivan Kisel  //
24 //                                                                           //
25 ///////////////////////////////////////////////////////////////////////////////
26
27 #if __GNUC__>= 3
28 using namespace std;
29 #endif
30
31 #include "AliHLTTPCCATrackerComponent.h"
32 #include "AliHLTTPCTransform.h"
33 #include "AliHLTTPCCATrackerFramework.h"
34 #include "AliHLTTPCCAOutTrack.h"
35 #include "AliHLTTPCCAParam.h"
36 #include "AliHLTTPCCATrackConvertor.h"
37 #include "AliHLTArray.h"
38
39 #include "AliHLTTPCSpacePointData.h"
40 #include "AliHLTTPCClusterDataFormat.h"
41 #include "AliHLTTPCCACompressedInputData.h"
42 #include "AliHLTTPCTransform.h"
43 #include "AliHLTTPCTrackSegmentData.h"
44 #include "AliHLTTPCTrackArray.h"
45 #include "AliHLTTPCTrackletDataFormat.h"
46 #include "AliHLTTPCDefinitions.h"
47 #include "AliExternalTrackParam.h"
48 #include "TStopwatch.h"
49 #include "TMath.h"
50 #include "AliCDBEntry.h"
51 #include "AliCDBManager.h"
52 #include "TObjString.h"
53 #include "TObjArray.h"
54 #include "AliHLTTPCCASliceOutput.h"
55 #include "AliHLTTPCCAClusterData.h"
56
57 const AliHLTComponentDataType AliHLTTPCCADefinitions::fgkTrackletsDataType = AliHLTComponentDataTypeInitializer( "CATRACKL", kAliHLTDataOriginTPC );
58
59 /** ROOT macro for the implementation of ROOT specific class methods */
60 ClassImp( AliHLTTPCCATrackerComponent )
61
62 AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent()
63     :
64     fTracker( NULL ),
65     fSolenoidBz( 0 ),
66     fMinNTrackClusters( 0 ),
67     fClusterZCut( 500. ),
68     fNeighboursSearchArea( 0 ), 
69     fClusterErrorCorrectionY(0), 
70     fClusterErrorCorrectionZ(0),
71     fFullTime( 0 ),
72     fRecoTime( 0 ),
73     fNEvents( 0 ),
74     fOutputTRAKSEGS( 0 ),
75     fAllowGPU( 0)
76 {
77   // see header file for class documentation
78   // or
79   // refer to README to build package
80   // or
81   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
82 }
83
84 AliHLTTPCCATrackerComponent::AliHLTTPCCATrackerComponent( const AliHLTTPCCATrackerComponent& )
85     :
86     AliHLTProcessor(),
87     fTracker( NULL ),
88     fSolenoidBz( 0 ),
89     fMinNTrackClusters( 30 ),
90     fClusterZCut( 500. ),
91     fNeighboursSearchArea(0),
92     fClusterErrorCorrectionY(0), 
93     fClusterErrorCorrectionZ(0),
94     fFullTime( 0 ),
95     fRecoTime( 0 ),
96     fNEvents( 0 ),
97     fOutputTRAKSEGS( 0 ),
98     fAllowGPU( 0)
99 {
100   // see header file for class documentation
101   HLTFatal( "copy constructor untested" );
102 }
103
104 AliHLTTPCCATrackerComponent& AliHLTTPCCATrackerComponent::operator=( const AliHLTTPCCATrackerComponent& )
105 {
106   // see header file for class documentation
107   HLTFatal( "assignment operator untested" );
108   return *this;
109 }
110
111 AliHLTTPCCATrackerComponent::~AliHLTTPCCATrackerComponent()
112 {
113   // see header file for class documentation
114   if (fTracker) delete fTracker;
115 }
116
117 //
118 // Public functions to implement AliHLTComponent's interface.
119 // These functions are required for the registration process
120 //
121
122 const char* AliHLTTPCCATrackerComponent::GetComponentID()
123 {
124   // see header file for class documentation
125   return "TPCCATracker";
126 }
127
128 void AliHLTTPCCATrackerComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list )
129 {
130   // see header file for class documentation
131   list.clear();
132   list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
133   list.push_back( AliHLTTPCCADefinitions::fgkCompressedInputDataType );
134 }
135
136 AliHLTComponentDataType AliHLTTPCCATrackerComponent::GetOutputDataType()
137 {
138   // see header file for class documentation
139   if ( fOutputTRAKSEGS ) return AliHLTTPCDefinitions::fgkTrackSegmentsDataType;
140   else return AliHLTTPCCADefinitions::fgkTrackletsDataType;
141 }
142
143 void AliHLTTPCCATrackerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
144 {
145   // define guess for the output data size
146   constBase = 200;       // minimum size
147   inputMultiplier = 3.; // size relative to input
148 }
149
150 AliHLTComponent* AliHLTTPCCATrackerComponent::Spawn()
151 {
152   // see header file for class documentation
153   return new AliHLTTPCCATrackerComponent;
154 }
155
156 void AliHLTTPCCATrackerComponent::SetDefaultConfiguration()
157 {
158   // Set default configuration for the CA tracker component
159   // Some parameters can be later overwritten from the OCDB
160
161   fSolenoidBz = -5.00668;
162   fMinNTrackClusters = 0;
163   fClusterZCut = 500.;
164   fNeighboursSearchArea = 0;
165   fClusterErrorCorrectionY = 0;
166   fClusterErrorCorrectionZ = 0;
167   fOutputTRAKSEGS = 0;
168   fFullTime = 0;
169   fRecoTime = 0;
170   fNEvents = 0;
171 }
172
173 int AliHLTTPCCATrackerComponent::ReadConfigurationString(  const char* arguments )
174 {
175   // Set configuration parameters for the CA tracker component from the string
176
177   int iResult = 0;
178   if ( !arguments ) return iResult;
179
180   TString allArgs = arguments;
181   TString argument;
182   int bMissingParam = 0;
183
184   TObjArray* pTokens = allArgs.Tokenize( " " );
185
186   int nArgs =  pTokens ? pTokens->GetEntries() : 0;
187
188   for ( int i = 0; i < nArgs; i++ ) {
189     argument = ( ( TObjString* )pTokens->At( i ) )->GetString();
190     if ( argument.IsNull() ) continue;
191
192     if ( argument.CompareTo( "-solenoidBz" ) == 0 ) {
193       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
194       HLTWarning("argument -solenoidBz is deprecated, magnetic field set up globally (%f)", GetBz());
195       continue;
196     }
197
198     if ( argument.CompareTo( "-minNClustersOnTrack" ) == 0 ) {
199       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
200       fMinNTrackClusters = ( ( TObjString* )pTokens->At( i ) )->GetString().Atoi();
201       HLTInfo( "minNClustersOnTrack set to: %d", fMinNTrackClusters );
202       continue;
203     }
204
205     if ( argument.CompareTo( "-clusterZCut" ) == 0 ) {
206       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
207       fClusterZCut = TMath::Abs( ( ( TObjString* )pTokens->At( i ) )->GetString().Atof() );
208       HLTInfo( "ClusterZCut set to: %f", fClusterZCut );
209       continue;
210     }
211  
212    if ( argument.CompareTo( "-neighboursSearchArea" ) == 0 ) {
213       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
214       fNeighboursSearchArea = TMath::Abs( ( ( TObjString* )pTokens->At( i ) )->GetString().Atof() );
215       HLTInfo( "NeighboursSearchArea set to: %f", fNeighboursSearchArea );
216       continue;
217     }
218
219    if ( argument.CompareTo( "-errorCorrectionY" ) == 0 ) {
220      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
221      fClusterErrorCorrectionY = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
222      HLTInfo( "Cluster Y error correction factor set to: %f", fClusterErrorCorrectionY );
223      continue;
224    }
225    
226    if ( argument.CompareTo( "-errorCorrectionZ" ) == 0 ) {
227      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
228      fClusterErrorCorrectionZ = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
229      HLTInfo( "Cluster Z error correction factor set to: %f", fClusterErrorCorrectionZ );
230      continue;
231    }
232
233    if ( argument.CompareTo( "-outputTRAKSEGS" ) == 0 ) {
234       fOutputTRAKSEGS = 1;
235       HLTInfo( "The special output type \"TRAKSEGS\" is set" );
236       continue;
237     }
238
239     if (argument.CompareTo( "-allowGPU" ) == 0) {
240       fAllowGPU = 1;
241       HLTImportant( "Will try to run tracker on GPU" );
242       continue;
243     }
244
245     HLTError( "Unknown option \"%s\"", argument.Data() );
246     iResult = -EINVAL;
247   }
248   delete pTokens;
249
250   if ( bMissingParam ) {
251     HLTError( "Specifier missed for parameter \"%s\"", argument.Data() );
252     iResult = -EINVAL;
253   }
254
255   return iResult;
256 }
257
258
259 int AliHLTTPCCATrackerComponent::ReadCDBEntry( const char* cdbEntry, const char* chainId )
260 {
261   // see header file for class documentation
262
263   const char* defaultNotify = "";
264
265   if ( !cdbEntry ) {
266     cdbEntry = "HLT/ConfigTPC/TPCCATracker";
267     defaultNotify = " (default)";
268     chainId = 0;
269   }
270
271   HLTInfo( "configure from entry \"%s\"%s, chain id %s", cdbEntry, defaultNotify, ( chainId != NULL && chainId[0] != 0 ) ? chainId : "<none>" );
272   AliCDBEntry *pEntry = AliCDBManager::Instance()->Get( cdbEntry );//,GetRunNo());
273
274   if ( !pEntry ) {
275     HLTError( "cannot fetch object \"%s\" from CDB", cdbEntry );
276     return -EINVAL;
277   }
278
279   TObjString* pString = dynamic_cast<TObjString*>( pEntry->GetObject() );
280
281   if ( !pString ) {
282     HLTError( "configuration object \"%s\" has wrong type, required TObjString", cdbEntry );
283     return -EINVAL;
284   }
285
286   HLTInfo( "received configuration object string: \"%s\"", pString->GetString().Data() );
287
288   return  ReadConfigurationString( pString->GetString().Data() );
289 }
290
291
292 int AliHLTTPCCATrackerComponent::Configure( const char* cdbEntry, const char* chainId, const char *commandLine )
293 {
294   // Configure the component
295   // There are few levels of configuration,
296   // parameters which are set on one step can be overwritten on the next step
297
298   //* read hard-coded values
299
300   SetDefaultConfiguration();
301
302   //* read the default CDB entry
303
304   int iResult1 = ReadCDBEntry( NULL, chainId );
305
306   //* read magnetic field
307
308   int iResult2 = 0; //ReadCDBEntry( kAliHLTCDBSolenoidBz, chainId );
309   fSolenoidBz = GetBz();
310
311   //* read the actual CDB entry if required
312
313   int iResult3 = ( cdbEntry ) ? ReadCDBEntry( cdbEntry, chainId ) : 0;
314
315   //* read extra parameters from input (if they are)
316
317   int iResult4 = 0;
318
319   if ( commandLine && commandLine[0] != '\0' ) {
320     HLTInfo( "received configuration string from HLT framework: \"%s\"", commandLine );
321     iResult4 = ReadConfigurationString( commandLine );
322   }
323
324   // Initialise the tracker here
325
326   return iResult1 ? iResult1 : ( iResult2 ? iResult2 : ( iResult3 ? iResult3 : iResult4 ) );
327 }
328
329
330
331 int AliHLTTPCCATrackerComponent::DoInit( int argc, const char** argv )
332 {
333   // Configure the CA tracker component
334
335   if ( fTracker ) return EINPROGRESS;
336
337
338   //fTracker = new AliHLTTPCCATrackerFramework();
339   //Do not initialize the TrackerFramework here since the CUDA framework is thread local and DoInit is called from different thread than DoEvent
340
341   TString arguments = "";
342   for ( int i = 0; i < argc; i++ ) {
343     if ( !arguments.IsNull() ) arguments += " ";
344     arguments += argv[i];
345   }
346
347   return Configure( NULL, NULL, arguments.Data() );
348 }
349
350
351 int AliHLTTPCCATrackerComponent::DoDeinit()
352 {
353   // see header file for class documentation
354   if (fTracker) delete fTracker;
355   fTracker = NULL;
356   return 0;
357 }
358
359
360
361 int AliHLTTPCCATrackerComponent::Reconfigure( const char* cdbEntry, const char* chainId )
362 {
363   // Reconfigure the component from OCDB .
364
365   return Configure( cdbEntry, chainId, NULL );
366 }
367
368
369 bool AliHLTTPCCATrackerComponent::CompareClusters( AliHLTTPCSpacePointData *a, AliHLTTPCSpacePointData *b )
370 {
371   //* Comparison function for sort clusters
372
373   if ( a->fPadRow < b->fPadRow ) return 1;
374   if ( a->fPadRow > b->fPadRow ) return 0;
375   return ( a->fZ < b->fZ );
376 }
377
378
379
380 int AliHLTTPCCATrackerComponent::DoEvent
381 (
382   const AliHLTComponentEventData& evtData,
383   const AliHLTComponentBlockData* blocks,
384   AliHLTComponentTriggerData& /*trigData*/,
385   AliHLTUInt8_t* outputPtr,
386   AliHLTUInt32_t& size,
387   vector<AliHLTComponentBlockData>& outputBlocks )
388 {
389   //* process event
390
391   AliHLTUInt32_t maxBufferSize = size;
392   size = 0; // output size
393
394   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) ) {
395     return 0;
396   }
397
398   TStopwatch timer;
399
400   // Event reconstruction in one TPC slice with CA Tracker
401
402   //Logging( kHLTLogWarning, "HLT::TPCCATracker::DoEvent", "DoEvent", "CA::DoEvent()" );
403   if ( evtData.fBlockCnt <= 0 ) {
404     HLTWarning( "no blocks in event" );
405     return 0;
406   }
407
408   const AliHLTComponentBlockData* iter = NULL;
409   unsigned long ndx;
410
411   // Determine the slice number
412
413   //Find min and max slice number with now slice missing in between (for default output)
414   int minslice = -1, maxslice = -1;
415   int slice = -1;
416   {
417     std::vector<int> slices;
418     std::vector<int>::iterator slIter;
419     std::vector<unsigned> sliceCnts;
420     std::vector<unsigned>::iterator slCntIter;
421
422     for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ) {
423       iter = blocks + ndx;
424       if ( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType
425            && iter->fDataType != AliHLTTPCCADefinitions::fgkCompressedInputDataType
426            ) continue;
427
428       slice = AliHLTTPCDefinitions::GetMinSliceNr( *iter );
429           if (slice < minslice || minslice == -1) minslice = slice;
430           if (slice > maxslice) maxslice = slice;
431
432       bool found = 0;
433       slCntIter = sliceCnts.begin();
434       for ( slIter = slices.begin(); slIter != slices.end(); slIter++, slCntIter++ ) {
435         if ( *slIter == slice ) {
436           found = kTRUE;
437           break;
438         }
439       }
440       if ( !found ) {
441         slices.push_back( slice );
442         sliceCnts.push_back( 1 );
443       } else (*slCntIter)++;
444     }
445
446           if ( slices.size() == 0 ) {
447                 HLTWarning( "no slices found in event" );
448                 return 0;
449           }
450
451
452     // Determine slice number to really use. (for obsolete output)
453     if ( slices.size() > 1 ) {
454                 Logging( fOutputTRAKSEGS ? kHLTLogError : kHLTLogDebug, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
455                "Multiple slice numbers found in event 0x%08lX (%lu). Determining maximum occuring slice number...",
456                evtData.fEventID, evtData.fEventID );
457       unsigned maxCntSlice = 0;
458       slCntIter = sliceCnts.begin();
459       for ( slIter = slices.begin(); slIter != slices.end(); slIter++, slCntIter++ ) {
460         Logging( fOutputTRAKSEGS ? kHLTLogError : kHLTLogDebug, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
461                  "Slice %lu found %lu times.", *slIter, *slCntIter );
462         if ( fOutputTRAKSEGS && maxCntSlice < *slCntIter ) {
463           maxCntSlice = *slCntIter;
464           slice = *slIter;
465         }
466       }
467       if (fOutputTRAKSEGS)
468       {
469             Logging( kHLTLogError, "HLT::TPCSliceTracker::DoEvent", "Multiple slices found in event",
470                "Using slice %lu.", slice );
471         }
472     } else if ( slices.size() > 0 ) {
473       slice = *( slices.begin() );
474     }
475
476
477           if (fOutputTRAKSEGS)
478           {
479                   minslice = maxslice = slice;
480           }
481           else
482           {
483                   for (int islice = minslice;islice <= maxslice;islice++)
484                   {
485                           bool found = false;
486                           for(slIter = slices.begin(); slIter != slices.end();slIter++)
487                           {
488                                   if (*slIter == islice)
489                                   {
490                                           found = true;
491                                           break;
492                                   }
493                           }
494                           if (!found)
495                           {
496                                   maxslice = islice - 1;
497                                   break;
498                           }
499                   }
500          }
501   }
502
503   if ( !fTracker ) fTracker = new AliHLTTPCCATrackerFramework(fAllowGPU);
504
505   int slicecount = maxslice + 1 - minslice;
506   if (slicecount > fTracker->MaxSliceCount())
507   {
508         maxslice = minslice + (slicecount = fTracker->MaxSliceCount()) - 1;
509   }
510   int nClustersTotalSum = 0;
511   AliHLTTPCCAClusterData* clusterData = new AliHLTTPCCAClusterData[slicecount];
512
513
514   // min and max patch numbers and row numbers
515   int* slicerow = new int[slicecount * 2];
516   int* sliceminPatch = new int[slicecount];
517   int* slicemaxPatch = new int[slicecount];
518   memset(slicerow, 0, slicecount * 2 * sizeof(int));
519   for (int i = 0;i < slicecount;i++)
520   {
521           sliceminPatch[i] = 100;
522           slicemaxPatch[i] = -1;
523   }
524
525   //Prepare everything for all slices
526
527   for (int islice = 0;islice < slicecount;islice++)
528   {
529           slice = minslice + islice;
530
531           // Initialize the tracker slice
532           {
533                 int iSec = slice;
534                 float inRmin = 83.65;
535                 //    float inRmax = 133.3;
536                 //    float outRmin = 133.5;
537                 float outRmax = 247.7;
538                 float plusZmin = 0.0529937;
539                 float plusZmax = 249.778;
540                 float minusZmin = -249.645;
541                 float minusZmax = -0.0799937;
542                 float dalpha = 0.349066;
543                 float alpha = 0.174533 + dalpha * iSec;
544
545                 bool zPlus = ( iSec < 18 );
546                 float zMin =  zPlus ? plusZmin : minusZmin;
547                 float zMax =  zPlus ? plusZmax : minusZmax;
548                 //TPCZmin = -249.645, ZMax = 249.778
549                 //    float rMin =  inRmin;
550                 //    float rMax =  outRmax;
551                 int nRows = AliHLTTPCTransform::GetNRows();
552
553                 float padPitch = 0.4;
554                 float sigmaZ = 0.228808;
555
556                 float *rowX = new float [nRows];
557                 for ( int irow = 0; irow < nRows; irow++ ) {
558                   rowX[irow] = AliHLTTPCTransform::Row2X( irow );
559                 }
560
561                 AliHLTTPCCAParam param;
562
563                 param.Initialize( iSec, nRows, rowX, alpha, dalpha,
564                                                   inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, fSolenoidBz );
565                 param.SetHitPickUpFactor( 2 );
566                 if( fNeighboursSearchArea>0 ) param.SetNeighboursSearchArea( fNeighboursSearchArea );
567                 if( fClusterErrorCorrectionY>1.e-4 ) param.SetClusterError2CorrectionY( fClusterErrorCorrectionY*fClusterErrorCorrectionY );
568                 if( fClusterErrorCorrectionZ>1.e-4 ) param.SetClusterError2CorrectionZ( fClusterErrorCorrectionZ*fClusterErrorCorrectionZ );
569                 param.Update();
570                 fTracker->InitializeSliceParam( slice, param );
571                 delete[] rowX;
572           }
573
574           // total n Hits
575           int nClustersTotal = 0;
576
577           // sort patches
578           std::vector<unsigned long> patchIndices;
579
580           for ( ndx = 0; ndx < evtData.fBlockCnt; ndx++ ) {
581                 iter = blocks + ndx;
582                 if ( slice != AliHLTTPCDefinitions::GetMinSliceNr( *iter ) ) continue;
583                 if ( iter->fDataType == AliHLTTPCDefinitions::fgkClustersDataType ){
584                   AliHLTTPCClusterData* inPtrSP = ( AliHLTTPCClusterData* )( iter->fPtr );
585                   nClustersTotal += inPtrSP->fSpacePointCnt;
586                 } 
587                 else if ( iter->fDataType == AliHLTTPCCADefinitions::fgkCompressedInputDataType){
588                   const AliHLTUInt8_t * inPtr =  (const AliHLTUInt8_t *)iter->fPtr;
589                   while( inPtr< ((const AliHLTUInt8_t *) iter->fPtr) + iter->fSize ){
590                     AliHLTTPCCACompressedClusterRow *row = (AliHLTTPCCACompressedClusterRow*)inPtr;
591                     nClustersTotal+= row->fNClusters;     
592                     inPtr = (const AliHLTUInt8_t *)(row->fClusters+row->fNClusters);
593                   }
594                 }
595                 else continue;
596
597                 int patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
598                 if ( sliceminPatch[islice] > patch ) {
599                   sliceminPatch[islice] = patch;
600                   slicerow[2 * islice + 0] = AliHLTTPCTransform::GetFirstRow( patch );
601                 }
602                 if ( slicemaxPatch[islice] < patch ) {
603                   slicemaxPatch[islice] = patch;
604                   slicerow[2 * islice + 1] = AliHLTTPCTransform::GetLastRow( patch );
605                 }
606                 std::vector<unsigned long>::iterator pIter = patchIndices.begin();
607                 while ( pIter != patchIndices.end() && AliHLTTPCDefinitions::GetMinPatchNr( blocks[*pIter] ) < patch ) {
608                   pIter++;
609                 }
610                 patchIndices.insert( pIter, ndx );
611           }
612
613
614           // pass event to CA Tracker
615
616
617           Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reading hits",
618                            "Total %d hits to read for slice %d", nClustersTotal, slice );
619
620
621           clusterData[islice].StartReading( slice, nClustersTotal );
622
623           for ( std::vector<unsigned long>::iterator pIter = patchIndices.begin(); pIter != patchIndices.end(); pIter++ ) {
624                 ndx = *pIter;
625                 iter = blocks + ndx;
626                 int patch = AliHLTTPCDefinitions::GetMinPatchNr( *iter );
627                 int nPatchClust = 0;
628                 if ( iter->fDataType == AliHLTTPCDefinitions::fgkClustersDataType ){
629                   AliHLTTPCClusterData* inPtrSP = ( AliHLTTPCClusterData* )( iter->fPtr );
630                   nPatchClust = inPtrSP->fSpacePointCnt;
631                   for ( unsigned int i = 0; i < inPtrSP->fSpacePointCnt; i++ ) {
632                     AliHLTTPCSpacePointData *c = &( inPtrSP->fSpacePoints[i] );
633                     if ( CAMath::Abs( c->fZ ) > fClusterZCut ) continue;
634                     if ( c->fPadRow > 159 ) {
635                       HLTError( "Wrong TPC cluster with row number %d received", c->fPadRow );
636                       continue;
637                     }
638                     clusterData[islice].ReadCluster( c->fID, c->fPadRow, c->fX, c->fY, c->fZ, c->fCharge );
639                   }           
640                 } 
641                 else if ( iter->fDataType == AliHLTTPCCADefinitions::fgkCompressedInputDataType){
642                   const AliHLTUInt8_t * inPtr = (const AliHLTUInt8_t *)iter->fPtr;
643                   nPatchClust=0;
644                   while( inPtr< ((const AliHLTUInt8_t *)iter->fPtr) + iter->fSize ){
645                     AliHLTTPCCACompressedClusterRow *row = (AliHLTTPCCACompressedClusterRow*)inPtr;
646                     UInt_t id = row->fSlicePatchRowID;
647                     UInt_t jslice = id>>10;       
648                     UInt_t jpatch = (id>>6) & 0x7;
649                     UInt_t jrow   =  id     & 0x3F;     
650                     jrow+= AliHLTTPCTransform::GetFirstRow( jpatch );
651                     Double_t rowX = AliHLTTPCTransform::Row2X( jrow );
652                     //cout<<"Read row: s "<<jslice<<" p "<<jpatch<<" r "<<jrow<<" x "<<row->fX<<" nclu "<<row->fNClusters<<" :"<<endl;
653                     if( jrow > 159 ) {
654                       HLTError( "Wrong TPC cluster with row number %d received", jrow );
655                       continue;
656                     }
657                     for ( unsigned int i = 0; i < row->fNClusters; i++ ) {
658                       AliHLTTPCCACompressedCluster *c = &( row->fClusters[i] );
659                       
660                       UInt_t ix0 = c->fP0 >>24;
661                       UInt_t ix1 = c->fP1 >>24;
662                       Double_t x = (ix1<<8) + ix0;
663                       Double_t y = c->fP0 & 0x00FFFFFF;
664                       Double_t z = c->fP1 & 0x00FFFFFF;
665                       x = (x - 32768.)*1.e-4 + rowX;
666                       y = (y - 8388608.)*1.e-4;
667                       z = (z - 8388608.)*1.e-4;
668                       
669                       UInt_t cluId = nPatchClust + ((jslice&0x7f)<<25)+((jpatch&0x7)<<22);
670                       //cout<<"clu "<<i<<": "<<x<<" "<<y<<" "<<z<<" "<<cluId<<endl;
671                       if ( CAMath::Abs( z ) <= fClusterZCut ){
672                         clusterData[islice].ReadCluster( cluId, jrow, x, y, z, 0 );
673                       }
674                       nPatchClust++;              
675                     }
676                     inPtr = (const AliHLTUInt8_t *)(row->fClusters+row->fNClusters);
677                   }
678                 }
679                 Logging( kHLTLogInfo, "HLT::TPCCATracker::DoEvent", "Reading hits",
680                          "Read %d hits for slice %d - patch %d", nPatchClust, slice, patch );
681           }
682
683           clusterData[islice].FinishReading();
684           nClustersTotalSum += nClustersTotal;
685   }
686
687   //Prepare Output
688   AliHLTTPCCASliceOutput::outputControlStruct outputControl;
689   //Set tracker output so tracker does not have to output both formats!
690   outputControl.fEndOfSpace = 0;
691   outputControl.fObsoleteOutput = fOutputTRAKSEGS;
692   outputControl.fDefaultOutput = !fOutputTRAKSEGS;
693
694   //For new output we can write directly to output buffer
695   outputControl.fOutputPtr = fOutputTRAKSEGS ? NULL : (char*) outputPtr;
696   outputControl.fOutputMaxSize = maxBufferSize;
697
698   AliHLTTPCCASliceOutput** sliceOutput = new AliHLTTPCCASliceOutput*[slicecount];
699   memset(sliceOutput, 0, slicecount * sizeof(AliHLTTPCCASliceOutput*));
700
701   // reconstruct the event
702   TStopwatch timerReco;
703   fTracker->SetOutputControl(&outputControl);
704   fTracker->ProcessSlices(minslice, slicecount, clusterData, sliceOutput);
705   timerReco.Stop();
706   
707   int ret = 0;
708   unsigned int mySize = 0;
709   int ntracks = 0;
710   int error = 0;
711
712   for (int islice = 0;islice < slicecount;islice++)
713   {
714     if( outputControl.fEndOfSpace ){
715       HLTWarning( "Output buffer size exceed, tracks are not stored" );
716       ret = -ENOSPC;
717       error = 1;
718       break;     
719     }
720           slice = minslice + islice;
721
722           if (sliceOutput[islice])
723           {
724                   // write reconstructed tracks
725
726                   if ( fOutputTRAKSEGS ) {
727
728                   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reconstruct",
729                                    "%d tracks found for slice %d", sliceOutput[islice]->NOutTracks(), slice );
730
731                         ntracks = sliceOutput[islice]->NOutTracks();
732
733                         AliHLTTPCTrackletData* outPtr = ( AliHLTTPCTrackletData* )( outputPtr );
734
735                         AliHLTTPCTrackSegmentData* currOutTracklet = outPtr->fTracklets;
736
737                         mySize =   ( ( AliHLTUInt8_t * )currOutTracklet ) -  ( ( AliHLTUInt8_t * )outputPtr );
738
739                         outPtr->fTrackletCnt = 0;
740
741                         for ( int itr = 0; itr < ntracks; itr++ ) {
742
743                           AliHLTTPCCAOutTrack &t = sliceOutput[islice]->OutTracks()[itr];
744
745                           //Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Wrtite output","track %d with %d hits", itr, t.NHits());
746
747                           if ( t.NHits() < fMinNTrackClusters ) continue;
748
749                           // calculate output track size
750
751                           unsigned int dSize = sizeof( AliHLTTPCTrackSegmentData ) + t.NHits() * sizeof( unsigned int );
752
753                           if ( mySize + dSize > maxBufferSize ) {
754                                 HLTWarning( "Output buffer size exceed (buffer size %d, current size %d), %d tracks are not stored", maxBufferSize, mySize, ntracks - itr + 1 );
755                                 ret = -ENOSPC;
756                                 error = 1;
757                                 break;
758                           }
759
760                           // convert CA track parameters to HLT Track Segment
761
762                           int iFirstRow = 1000;
763                           int iLastRow = -1;
764                           int iFirstHit = sliceOutput[islice]->OutTrackHit(t.FirstHitRef());
765                           int iLastHit = iFirstHit;
766                           for ( int ih = 0; ih < t.NHits(); ih++ ) {
767                                 int hitID = sliceOutput[islice]->OutTrackHit(t.FirstHitRef() + ih);
768                                 int iRow = clusterData[islice].RowNumber( hitID );
769                                 if ( iRow < iFirstRow ) {  iFirstRow = iRow; iFirstHit = hitID; }
770                                 if ( iRow > iLastRow ) { iLastRow = iRow; iLastHit = hitID; }
771                           }
772
773                           AliHLTTPCCATrackParam par = t.StartPoint();
774
775                           par.TransportToX( clusterData[islice].X( iFirstHit ), .99 );
776
777                           AliExternalTrackParam tp;
778                           AliHLTTPCCATrackConvertor::GetExtParam( par, tp, 0 );
779
780                           currOutTracklet->fX = tp.GetX();
781                           currOutTracklet->fY = tp.GetY();
782                           currOutTracklet->fZ = tp.GetZ();
783                           currOutTracklet->fCharge = ( int ) tp.GetSign();
784                           currOutTracklet->fPt = TMath::Abs( tp.GetSignedPt() );
785                           float snp =  tp.GetSnp() ;
786                           if ( snp > .999 ) snp = .999;
787                           if ( snp < -.999 ) snp = -.999;
788                           currOutTracklet->fPsi = TMath::ASin( snp );
789                           currOutTracklet->fTgl = tp.GetTgl();
790
791                           currOutTracklet->fY0err = tp.GetSigmaY2();
792                           currOutTracklet->fZ0err = tp.GetSigmaZ2();
793                           float h = -currOutTracklet->fPt * currOutTracklet->fPt;
794                           currOutTracklet->fPterr = h * h * tp.GetSigma1Pt2();
795                           h = 1. / TMath::Sqrt((1.-snp)*(1.+snp));
796                           currOutTracklet->fPsierr = h * h * tp.GetSigmaSnp2();
797                           currOutTracklet->fTglerr = tp.GetSigmaTgl2();
798
799                           if ( par.TransportToX( clusterData[islice].X( iLastHit ), .99 ) ) {
800                                 currOutTracklet->fLastX = par.GetX();
801                                 currOutTracklet->fLastY = par.GetY();
802                                 currOutTracklet->fLastZ = par.GetZ();
803                           } else {
804                                 currOutTracklet->fLastX = clusterData[islice].X( iLastHit );
805                                 currOutTracklet->fLastY = clusterData[islice].Y( iLastHit );
806                                 currOutTracklet->fLastZ = clusterData[islice].Z( iLastHit );
807                           }
808                           //if( currOutTracklet->fLastX<10. ) {
809                           //HLTError("CA last point: hitxyz=%f,%f,%f, track=%f,%f,%f, tracklet=%f,%f,%f, nhits=%d",clusterData[islice].X( iLastHit ),clusterData[islice].Y( iLastHit],clusterData[islice].Z( iLastHit],
810                           //par.GetX(), par.GetY(),par.GetZ(),currOutTracklet->fLastX,currOutTracklet->fLastY ,currOutTracklet->fLastZ, t.NHits());
811                           //}
812                 #ifdef INCLUDE_TPC_HOUGH
813                 #ifdef ROWHOUGHPARAMS
814                           currOutTracklet->fTrackID = 0;
815                           currOutTracklet->fRowRange1 = clusterData[islice].RowNumber( iFirstHit );
816                           currOutTracklet->fRowRange2 = clusterData[islice].RowNumber( iLastHit );
817                           currOutTracklet->fSector = slice;
818                           currOutTracklet->fPID = 211;
819                 #endif
820                 #endif // INCLUDE_TPC_HOUGH
821
822
823                           currOutTracklet->fNPoints = t.NHits();
824
825                           for ( int i = 0; i < t.NHits(); i++ ) {
826                                 currOutTracklet->fPointIDs[i] = clusterData[islice].Id( sliceOutput[islice]->OutTrackHit(t.FirstHitRef()+i) );
827                           }
828
829                           currOutTracklet = ( AliHLTTPCTrackSegmentData* )( ( Byte_t * )currOutTracklet + dSize );
830                           mySize += dSize;
831                           outPtr->fTrackletCnt++;
832                         }
833
834                   } else { // default output type
835
836                   Logging( kHLTLogDebug, "HLT::TPCCATracker::DoEvent", "Reconstruct",
837                                    "%d tracks found for slice %d", sliceOutput[islice]->NTracks(), slice );
838
839                           mySize += sliceOutput[islice]->OutputMemorySize();
840                           ntracks += sliceOutput[islice]->NTracks();
841                   }
842           }
843           else
844           {
845                   HLTWarning( "Output buffer size exceed (buffer size %d, current size %d), tracks are not stored", maxBufferSize, mySize );
846                   mySize = 0;
847                   ret = -ENOSPC;
848                   ntracks = 0;
849                   error = 1;
850                   break;
851           }
852   }
853
854   size = 0;
855   if (error == 0)
856   {
857           for (int islice = 0;islice < slicecount;islice++)
858           {
859                   slice = minslice + islice;
860                   if (!fOutputTRAKSEGS) mySize = sliceOutput[islice]->OutputMemorySize();
861                   if (mySize > 0)
862                   {
863                         AliHLTComponentBlockData bd;
864                         FillBlockData( bd );
865                         bd.fOffset = fOutputTRAKSEGS ? 0 : ((char*) sliceOutput[islice] - (char*) outputPtr);
866                         bd.fSize = mySize;
867                         bd.fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( slice, slice, sliceminPatch[islice], slicemaxPatch[islice] );
868                         bd.fDataType = GetOutputDataType();
869                         outputBlocks.push_back( bd );
870                         size += mySize;
871                   }
872           }
873   }
874
875   //No longer needed
876
877   delete[] clusterData;
878   delete[] sliceOutput;
879
880   timer.Stop();
881
882   fFullTime += timer.RealTime();
883   fRecoTime += timerReco.RealTime();
884   fNEvents++;
885
886   // Set log level to "Warning" for on-line system monitoring
887   int hz = ( int ) ( fFullTime > 1.e-10 ? fNEvents / fFullTime : 100000 );
888   int hz1 = ( int ) ( fRecoTime > 1.e-10 ? fNEvents / fRecoTime : 100000 );
889   //Min and Max Patch are taken for first slice processed...
890   HLTInfo( "CATracker slices %d-%d: output %d tracks;  input %d clusters, patches %d..%d, rows %d..%d; time: full %d / reco %d Hz",
891            minslice, maxslice, ntracks, nClustersTotalSum, sliceminPatch[0], slicemaxPatch[0], slicerow[0], slicerow[1], hz, hz1 );
892
893   //No longer needed
894
895   delete[] slicerow;
896   delete[] sliceminPatch;
897   delete[] slicemaxPatch;
898
899   return ret;
900 }
901
902