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