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