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