MC labels added for ITS tracks ?\127 ?\127?\127
[u/mrichter/AliRoot.git] / HLT / ITS / AliHLTITSVertexerSPDComponent.cxx
1 // $Id: AliHLTITSVertexerSPDComponent.cxx 32659 2009-06-02 16:08:40Z sgorbuno $
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 ///  @file   AliHLTITSVertexerSPDComponent.cxx
21 ///  @author Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>
22 ///  @date   June 2009
23 ///  @brief  An ITS tracker processing component for the HLT
24
25
26 /////////////////////////////////////////////////////
27 //                                                 //
28 // a ITS tracker processing component for the HLT  //
29 //                                                 //
30 /////////////////////////////////////////////////////
31
32 #if __GNUC__>= 3
33 using namespace std;
34 #endif
35
36 #include "AliHLTITSVertexerSPDComponent.h"
37 #include "AliHLTArray.h"
38 #include "AliExternalTrackParam.h"
39 #include "TStopwatch.h"
40 #include "TMath.h"
41 #include "AliCDBEntry.h"
42 #include "AliCDBManager.h"
43 #include "AliGeomManager.h"
44 #include "TObjString.h"
45 #include "TObjArray.h"
46 #include "AliITStrackerHLT.h"
47 #include "AliHLTITSSpacePointData.h"
48 #include "AliHLTITSClusterDataFormat.h"
49 #include "AliHLTDataTypes.h"
50 #include "AliHLTExternalTrackParam.h"
51 #include "AliHLTGlobalBarrelTrack.h"
52 #include "AliGeomManager.h"
53 #include "AliESDVertex.h"
54
55 /** ROOT macro for the implementation of ROOT specific class methods */
56 ClassImp( AliHLTITSVertexerSPDComponent )
57 AliHLTITSVertexerSPDComponent::AliHLTITSVertexerSPDComponent()
58     :
59     fZRange(40),
60     fZBinSize(1),
61     fFullTime( 0 ),
62     fRecoTime( 0 ),
63     fNEvents( 0 ),
64     fSumW(0),
65     fSumN(0),
66     fZMin(0),
67     fNZBins(0)
68 {
69   // see header file for class documentation
70   // or
71   // refer to README to build package
72   // or
73   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
74
75   for( int i=0; i<9; i++ ) fSum[i] = 0;
76
77   fRunVtx[0] = 0;
78   fRunVtx[1] = 0;
79   fRunVtx[2] = 0;
80 }
81
82 AliHLTITSVertexerSPDComponent::AliHLTITSVertexerSPDComponent( const AliHLTITSVertexerSPDComponent& )
83     :
84     AliHLTProcessor(),
85     fZRange(40),
86     fZBinSize(1),
87     fFullTime( 0 ),
88     fRecoTime( 0 ),
89     fNEvents( 0 ),
90     fSumW(0),
91     fSumN(0),
92     fZMin(0),
93     fNZBins(0)
94 {
95   // see header file for class documentation
96   HLTFatal( "copy constructor untested" );
97 }
98
99 AliHLTITSVertexerSPDComponent& AliHLTITSVertexerSPDComponent::operator=( const AliHLTITSVertexerSPDComponent& )
100 {
101   // see header file for class documentation
102   HLTFatal( "assignment operator untested" );
103   return *this;
104 }
105
106 AliHLTITSVertexerSPDComponent::~AliHLTITSVertexerSPDComponent()
107 {
108   // see header file for class documentation  
109   for( int i=0; i<9; i++ ){
110     delete[] fSum[i];
111     fSum[i] = 0;
112   }
113   delete[] fSumW;
114   delete[] fSumN;
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* AliHLTITSVertexerSPDComponent::GetComponentID()
123 {
124   // see header file for class documentation
125   return "ITSVertexerSPD";
126 }
127
128 void AliHLTITSVertexerSPDComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list )
129 {
130   // see header file for class documentation
131   list.clear();
132   list.push_back( kAliHLTDataTypeClusters|kAliHLTDataOriginITSSPD );
133 }
134
135 AliHLTComponentDataType AliHLTITSVertexerSPDComponent::GetOutputDataType()
136 {
137   // see header file for class documentation  
138   return kAliHLTMultipleDataType;
139 }
140
141 int AliHLTITSVertexerSPDComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
142
143 {
144   // see header file for class documentation
145   tgtList.clear();
146   tgtList.push_back( kAliHLTDataTypeESDVertex|kAliHLTDataOriginITS);
147   return tgtList.size();
148 }
149
150 void AliHLTITSVertexerSPDComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
151 {
152   // define guess for the output data size
153   constBase = 10000;       // minimum size
154   inputMultiplier = 0.5; // size relative to input
155 }
156
157 AliHLTComponent* AliHLTITSVertexerSPDComponent::Spawn()
158 {
159   // see header file for class documentation
160   return new AliHLTITSVertexerSPDComponent;
161 }
162
163 void AliHLTITSVertexerSPDComponent::SetDefaultConfiguration()
164 {
165   // Set default configuration for the CA tracker component
166   // Some parameters can be later overwritten from the OCDB
167
168   fRunVtx[0] = 0;
169   fRunVtx[1] = 0;
170   fRunVtx[2] = 0;
171   fZRange = 40;
172   fZBinSize = 1;
173   fFullTime = 0;
174   fRecoTime = 0;
175   fNEvents = 0;
176 }
177
178 int AliHLTITSVertexerSPDComponent::ReadConfigurationString(  const char* arguments )
179 {
180   // Set configuration parameters for the CA tracker component from the string
181
182   int iResult = 0;
183   if ( !arguments ) return iResult;
184
185   TString allArgs = arguments;
186   TString argument;
187   int bMissingParam = 0;
188
189   TObjArray* pTokens = allArgs.Tokenize( " " );
190
191   int nArgs =  pTokens ? pTokens->GetEntries() : 0;
192
193   for ( int i = 0; i < nArgs; i++ ) {
194     argument = ( ( TObjString* )pTokens->At( i ) )->GetString();
195     if ( argument.IsNull() ) continue;
196
197     if ( argument.CompareTo( "-runVertex" ) == 0 ) {
198       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
199       fRunVtx[0] = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
200       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
201       fRunVtx[1] = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
202       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
203       fRunVtx[2] = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
204       HLTInfo( "Default run vertex is set to (%f,%f,%f)",fRunVtx[0],
205                fRunVtx[1],fRunVtx[2] );
206       continue;
207     }
208
209     if ( argument.CompareTo( "-zRange" ) == 0 ) {
210       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
211       fZRange = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
212       HLTInfo( "Z range for the vertex search is set to +-%f cm", fZRange );
213       continue;
214     }
215
216     if ( argument.CompareTo( "-zBinSize" ) == 0 ) {
217       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
218       fZBinSize = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
219       HLTInfo( "Size of the Z bin for the vertex search is set to %f cm", fZBinSize );
220       continue;
221     }
222
223     HLTError( "Unknown option \"%s\"", argument.Data() );
224     iResult = -EINVAL;
225   }
226   delete pTokens;
227
228   if ( bMissingParam ) {
229     HLTError( "Specifier missed for parameter \"%s\"", argument.Data() );
230     iResult = -EINVAL;
231   }
232
233   return iResult;
234 }
235
236
237 int AliHLTITSVertexerSPDComponent::ReadCDBEntry( const char* cdbEntry, const char* chainId )
238 {
239   // see header file for class documentation
240
241   const char* defaultNotify = "";
242
243   if ( !cdbEntry ) {
244     return 0;// need to add the HLT/ConfigITS/ITSTracker directory to cdb SG!!!
245     cdbEntry = "HLT/ConfigITS/ITSTracker";
246     defaultNotify = " (default)";
247     chainId = 0;
248   }
249
250   HLTInfo( "configure from entry \"%s\"%s, chain id %s", cdbEntry, defaultNotify, ( chainId != NULL && chainId[0] != 0 ) ? chainId : "<none>" );
251   AliCDBEntry *pEntry = AliCDBManager::Instance()->Get( cdbEntry );//,GetRunNo());
252
253   if ( !pEntry ) {
254     HLTError( "cannot fetch object \"%s\" from CDB", cdbEntry );
255     return -EINVAL;
256   }
257
258   TObjString* pString = dynamic_cast<TObjString*>( pEntry->GetObject() );
259
260   if ( !pString ) {
261     HLTError( "configuration object \"%s\" has wrong type, required TObjString", cdbEntry );
262     return -EINVAL;
263   }
264
265   HLTInfo( "received configuration object string: \"%s\"", pString->GetString().Data() );
266
267   return  ReadConfigurationString( pString->GetString().Data() );
268 }
269
270
271 int AliHLTITSVertexerSPDComponent::Configure( const char* cdbEntry, const char* chainId, const char *commandLine )
272 {
273   // Configure the component
274   // There are few levels of configuration,
275   // parameters which are set on one step can be overwritten on the next step
276
277   //* read hard-coded values
278
279   SetDefaultConfiguration();
280
281   //* read the default CDB entry
282
283   int iResult1 = ReadCDBEntry( NULL, chainId );
284
285   //* read the actual CDB entry if required
286
287   int iResult2 = ( cdbEntry ) ? ReadCDBEntry( cdbEntry, chainId ) : 0;
288   
289   //* read extra parameters from input (if they are)
290
291   int iResult3 = 0;
292
293   if ( commandLine && commandLine[0] != '\0' ) {
294     HLTInfo( "received configuration string from HLT framework: \"%s\"", commandLine );
295     iResult3 = ReadConfigurationString( commandLine );
296   }
297
298   // Initialise the tracker here
299
300   return iResult1 ? iResult1 : ( iResult2 ? iResult2 : iResult3 );
301 }
302
303
304
305 int AliHLTITSVertexerSPDComponent::DoInit( int argc, const char** argv )
306 {
307   // Configure the ITS tracker component
308   
309   SetDefaultConfiguration();
310
311   if(AliGeomManager::GetGeometry()==NULL){
312     AliGeomManager::LoadGeometry("");
313   }
314   AliGeomManager::ApplyAlignObjsFromCDB("ITS");
315   
316   TString arguments = "";
317   for ( int i = 0; i < argc; i++ ) {
318     if ( !arguments.IsNull() ) arguments += " ";
319     arguments += argv[i];
320   }
321
322   int ret = Configure( NULL, NULL, arguments.Data() );
323
324   for( int i=0; i<9; i++ ) delete[] fSum[i];
325   delete[] fSumW;
326   delete[] fSumN;
327
328   fZMin = -fZRange;
329   fNZBins = ( (int) (2*fZRange/fZBinSize ))+1;
330
331   for( int i=0; i<9; i++ ) fSum[i] = new double [fNZBins];
332
333   fSumW = new double [fNZBins];
334   fSumN = new int [fNZBins];
335
336   return ret;
337 }
338
339
340 int AliHLTITSVertexerSPDComponent::DoDeinit()
341 {
342   // see header file for class documentation
343
344   for( int i=0; i<9; i++ ){    
345     delete[] fSum[i];
346     fSum[i] = 0;
347   }
348   delete[] fSumW;
349   delete[] fSumN;
350   fSumW = 0;
351   fSumN = 0;
352   
353   SetDefaultConfiguration();
354
355   return 0;
356 }
357
358
359
360 int AliHLTITSVertexerSPDComponent::Reconfigure( const char* cdbEntry, const char* chainId )
361 {
362   // Reconfigure the component from OCDB .
363
364   return Configure( cdbEntry, chainId, NULL );
365 }
366
367
368 int AliHLTITSVertexerSPDComponent::DoEvent
369 (
370   const AliHLTComponentEventData& evtData,
371   const AliHLTComponentBlockData* blocks,
372   AliHLTComponentTriggerData& /*trigData*/,
373   AliHLTUInt8_t* /*outputPtr*/,
374   AliHLTUInt32_t& size,
375   vector<AliHLTComponentBlockData>& /*outputBlocks*/ )
376 {
377   //* process event
378
379   //AliHLTUInt32_t maxBufferSize = size;
380   size = 0; // output size
381
382   if (!IsDataEvent()) return 0;
383
384   if ( evtData.fBlockCnt <= 0 ) {
385     HLTWarning( "no blocks in event" );
386     return 0;
387   }
388
389   TStopwatch timer;
390
391   // Event reconstruction in ITS
392
393   int iResult=0;
394
395   int nBlocks = evtData.fBlockCnt;
396   int nClustersTotal = 0;
397  
398
399   const int kNPhiBins = 20;
400   const double kDPhi = TMath::TwoPi() / kNPhiBins;
401   double vtxX = fRunVtx[0], vtxY = fRunVtx[1], vtxZ = fRunVtx[2];
402
403   std::vector<AliHLTITSVertexerSPDComponent::AliHLTITSVZCluster> clusters[2][ kNPhiBins ];
404
405   for (int ndx=0; ndx<nBlocks && iResult>=0; ndx++) {
406
407     const AliHLTComponentBlockData* iter = blocks+ndx;
408  
409     // Read ITS SPD clusters
410
411     if ( iter->fDataType == (kAliHLTDataTypeClusters|kAliHLTDataOriginITSSPD) ){
412
413       AliHLTITSClusterData *inPtr=reinterpret_cast<AliHLTITSClusterData*>( iter->fPtr );
414       int nClusters = inPtr->fSpacePointCnt;
415       nClustersTotal+=nClusters;
416       for( int icl=0; icl<nClusters; icl++ ){   
417         AliHLTITSSpacePointData &d = inPtr->fSpacePoints[icl];
418         if( d.fLayer>1 ) continue;// SPD only;
419         if( d.fLayer<0 ) continue;// SPD only;
420         Int_t lab[4] = { d.fTracks[0], d.fTracks[1], d.fTracks[2], d.fIndex };
421         Int_t info[3] = { d.fNy, d.fNz, d.fLayer };
422         Float_t hit[6] = { d.fY, d.fZ, d.fSigmaY2, d.fSigmaZ2, d.fQ, d.fSigmaYZ };
423         if( d.fLayer==4 ) hit[5] = -hit[5];
424
425
426         AliITSRecPoint p( lab, hit, info );
427
428         if (!p.Misalign()){
429           HLTWarning("Can not misalign an SPD cluster");
430         }
431         Float_t xyz[3];      
432         p.GetGlobalXYZ(xyz);      
433
434         // get phi bin
435         double phi = TMath::ATan2( xyz[1]-vtxY, xyz[0]-vtxX );
436         if( phi<0 ) phi+=TMath::TwoPi();
437         int iphi = (int ) phi/kDPhi;
438         if( iphi<0 ) iphi = 0;
439         if( iphi>=kNPhiBins ) iphi = kNPhiBins-1;
440         AliHLTITSVZCluster c;
441         c.fX = xyz[0];
442         c.fY = xyz[1];
443         c.fZ = xyz[2];
444         clusters[d.fLayer][iphi].push_back( c );        
445      }   
446     }    
447   }// end read input blocks
448   
449
450   // Reconstruct the vertex
451
452   TStopwatch timerReco;
453     
454   double zScale = 1./fZBinSize;
455
456   // fit few times decreasing the radius cut
457
458   double vtxDR2[5] = { 1.*1., .5*.5, .4*.4, .4*.4, .2*.2};
459   double vtxDZ [5] = { 1., .5, .3, .3, .3};
460   double  maxW = 0;
461   int bestBin = -1;
462   
463   for( int iter = 0; iter<3; iter++ ){
464
465     bool doZBins = (iter<2);
466     int nSearchBins = doZBins ?fNZBins :1;
467     {
468       for( int i=0; i<nSearchBins; i++ ){
469         fSum[0][i] = 0;
470         fSum[1][i] = 0;
471         fSum[2][i] = 0;
472         fSum[3][i] = 0;
473         fSum[4][i] = 0;
474         fSum[5][i] = 0;
475         fSum[6][i] = 0;
476         fSum[7][i] = 0;
477         fSum[8][i] = 0;
478         fSumW[i] = 0;
479         fSumN[i] = 0;
480       }
481     }
482
483     for( int iPhi=0; iPhi<kNPhiBins; iPhi++ ){
484       int nCluUp = clusters[1][iPhi].size();    
485       int nCluDn = clusters[0][iPhi].size();    
486       for( int icUp=0; icUp<nCluUp; icUp++ ){
487         AliHLTITSVZCluster &cu = clusters[1][iPhi][icUp];
488         double x0 = cu.fX - vtxX;
489         double y0 = cu.fY - vtxY;
490         double z0 = cu.fZ;// - vtxZ;
491         double bestR2 = 1.e10;
492         int bestDn=-1, bestBinDn =0;
493         double bestP[3]={0,0,0};
494         for( int icDn=0; icDn<nCluDn; icDn++ ){
495           AliHLTITSVZCluster &cd = clusters[0][iPhi][icDn];
496           double x1 = cd.fX - vtxX;
497           double y1 = cd.fY - vtxY;
498           double z1 = cd.fZ;// - vtxZ;
499           double dx = x1 - x0;
500           double dy = y1 - y0;
501           double l2 = 1./(dx*dx + dy*dy);
502           double a = x1*y0 - x0*y1;
503           double r2 = a*a*l2; 
504           if( r2 > vtxDR2[iter] ) continue;
505           if( r2 > bestR2 ) continue;
506           //double xv = -dy*a*l2;
507           //double yv =  dx*a*l2;
508           double zv = ( (x1*z0-x0*z1)*dx + (y1*z0-y0*z1)*dy )*l2;
509           int zbin;
510           if( doZBins ){
511             zbin = (int)((zv - fZMin)*zScale);
512             if( zbin<0 || zbin>=fNZBins ) continue;       
513           } else {
514             zbin = 0;
515             if( TMath::Abs( zv - vtxZ ) > vtxDZ[ iter ] ) continue;
516           }
517           bestR2 = r2;
518           bestDn = icDn;
519           bestP[0] = x1;
520           bestP[1] = y1;
521           bestP[2] = z1;
522           bestBinDn = zbin;       
523         }
524         if( bestDn < 0 ) continue;
525
526         double dx = bestP[0] - x0;
527         double dy = bestP[1] - y0;
528         double dz = bestP[2] - z0;        
529         double w = 1./(1.+bestR2);
530           
531         // Equations:
532         //
533         // A_i*x + B_i*y + C_i*z = D_i
534         //
535         // Sum[0] = sum A_i*A_i
536         // Sum[1] = sum B_i*A_i
537         // Sum[2] = sum B_i*B_i
538         // Sum[3] = sum C_i*A_i
539         // Sum[4] = sum C_i*B_i
540         // Sum[5] = sum C_i*C_i   
541         // Sum[6] = sum A_i*D_i
542         // Sum[7] = sum B_i*D_i
543         // Sum[8] = sum C_i*D_i
544         
545         double n = w / TMath::Sqrt(dx*dx + dy*dy + dz*dz);        
546         dy*=n;
547         dx*=n;
548         dz*=n;
549         double a, b, c, d;
550         a = dy; b = -dx; c = 0; d = dy*x0 - dx*y0;
551         {
552           fSum[0][bestBinDn]+= a*a;
553           fSum[1][bestBinDn]+= b*a;
554           fSum[2][bestBinDn]+= b*b;
555           //fSum[3][bestBinDn]+= c*a;
556           //fSum[4][bestBinDn]+= c*b;
557           //fSum[5][bestBinDn]+= c*c;
558           fSum[6][bestBinDn]+= a*d;
559           fSum[7][bestBinDn]+= b*d;
560           //fSum[8][bestBinDn]+= c*d;
561         }
562         a = dz; b = 0; c = -dx; d = dz*x0 - dx*z0;
563         {
564           fSum[0][bestBinDn]+= a*a;
565           //fSum[1][bestBinDn]+= b*a;
566           //fSum[2][bestBinDn]+= b*b;
567           fSum[3][bestBinDn]+= c*a;
568           //fSum[4][bestBinDn]+= c*b;
569           fSum[5][bestBinDn]+= c*c;
570           fSum[6][bestBinDn]+= a*d;
571           //fSum[7][bestBinDn]+= b*d;
572           fSum[8][bestBinDn]+= c*d;
573         }
574         
575         a = 0; b = dz; c = -dy; d = dz*y0 - dy*z0;
576         {
577           //fSum[0][bestBinDn]+= a*a;
578           //fSum[1][bestBinDn]+= b*a;
579           fSum[2][bestBinDn]+= b*b;
580           //fSum[3][bestBinDn]+= c*a;
581           fSum[4][bestBinDn]+= c*b;
582           fSum[5][bestBinDn]+= c*c;
583           //fSum[6][bestBinDn]+= a*d;
584           fSum[7][bestBinDn]+= b*d;
585           fSum[8][bestBinDn]+= c*d;
586         }
587
588         a = dz; b = 0; c = -dx; d = dz*x0 - dx*z0;
589         {
590           fSum[0][bestBinDn]+= a*a;
591           //fSum[1][bestBinDn]+= b*a;
592           //fSum[2][bestBinDn]+= b*b;
593           fSum[3][bestBinDn]+= c*a;
594           //fSum[4][bestBinDn]+= c*b;
595           fSum[5][bestBinDn]+= c*c;
596           fSum[6][bestBinDn]+= a*d;
597           //fSum[7][bestBinDn]+= b*d;
598           fSum[8][bestBinDn]+= c*d;
599         }
600         
601         fSumW[bestBinDn]+=w;
602         fSumN[bestBinDn]+=1;
603       }    
604     }
605     
606     maxW = 0;  
607     bestBin = -1;
608     for( int i=0; i<nSearchBins; i++ ){
609       if( fSumN[i]>maxW ){
610         bestBin = i;
611         maxW = fSumN[i];
612       }
613     }
614     if( bestBin<0 || fSumN[bestBin] < 3 ){
615       bestBin = -1;
616       break;
617     }
618     
619     // calculate the vertex position in the best bin
620     Double_t w[6] = { fSum[0][bestBin],
621                       fSum[1][bestBin], fSum[2][bestBin],
622                       fSum[3][bestBin], fSum[4][bestBin], fSum[5][bestBin] };
623     Double_t wI[6];
624     {
625       wI[0] = w[2]*w[5] - w[4]*w[4];
626       wI[1] = w[3]*w[4] - w[1]*w[5];
627       wI[2] = w[0]*w[5] - w[3]*w[3];
628       wI[3] = w[1]*w[4] - w[2]*w[3];
629       wI[4] = w[1]*w[3] - w[0]*w[4];
630       wI[5] = w[0]*w[2] - w[1]*w[1];     
631       
632       Double_t s = ( w[0]*wI[0] + w[1]*wI[1] + w[3]*wI[3] );
633       
634       if( s<1.e-10 ){
635         bestBin = -1;
636         break;
637       }
638       s = 1./s;   
639       wI[0]*=s;
640       wI[1]*=s;
641       wI[2]*=s;
642       wI[3]*=s;
643       wI[4]*=s;
644       wI[5]*=s;
645     }    
646
647     vtxX += wI[0]*fSum[6][bestBin] + wI[1]*fSum[7][bestBin] + wI[3]*fSum[8][bestBin];
648     vtxY += wI[1]*fSum[6][bestBin] + wI[2]*fSum[7][bestBin] + wI[4]*fSum[8][bestBin];
649     vtxZ = wI[3]*fSum[6][bestBin] + wI[4]*fSum[7][bestBin] + wI[5]*fSum[8][bestBin];
650     //cout<<"SG: "<<iter<<": "<<vtxX<<" "<<vtxY<<" "<<vtxZ<<endl;
651   }
652
653   timerReco.Stop();
654   
655   // Fill output 
656
657   if( bestBin>=0 ){
658     double pos[3] = {vtxX, vtxY, vtxZ};
659     double s = 400.E-4;
660     double cov[6] = {s*s,0,s*s,0,0,s*s};
661     AliESDVertex v(pos, cov, 0, fSumN[bestBin]);
662     PushBack( (TObject*) &v, kAliHLTDataTypeESDVertex|kAliHLTDataOriginITS,0 );
663
664     //cout<<"ITSVertexerSPD: vtx found "<<vtxX<<" "<<vtxY<<" "<<vtxZ<<endl;
665   }
666
667   timer.Stop();
668   fFullTime += timer.RealTime();
669   fRecoTime += timerReco.RealTime();
670   fNEvents++;
671
672   // Set log level to "Warning" for on-line system monitoring
673   int hz = ( int ) ( fFullTime > 1.e-10 ? fNEvents / fFullTime : 100000 );
674   int hz1 = ( int ) ( fRecoTime > 1.e-10 ? fNEvents / fRecoTime : 100000 );
675   HLTInfo( "ITS Z Vertexer: input %d clusters; time: full %d / reco %d Hz",
676            nClustersTotal, hz, hz1 );
677
678   return iResult;
679 }