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 }