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 }