HLT/global/AliHLTGlobalVertexerComponent.cxx

   1 //**************************************************************************
   2 //* This file is property of and copyright by the ALICE HLT Project        *
   3 //* ALICE Experiment at CERN, All rights reserved.                         *
   4 //*                                                                        *
   5 //* Primary Authors: S.Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>    *
   6 //*                  for The ALICE HLT Project.                            *
   7 //*                                                                        *
   8 //* Permission to use, copy, modify and distribute this software and its   *
   9 //* documentation strictly for non-commercial purposes is hereby granted   *
  10 //* without fee, provided that the above copyright notice appears in all   *
  11 //* copies and that both the copyright notice and this permission notice   *
  12 //* appear in the supporting documentation. The authors make no claims     *
  13 //* about the suitability of this software for any purpose. It is          *
  14 //* provided "as is" without express or implied warranty.                  *
  15 //**************************************************************************
  16
  17 /** @file   AliHLTGlobalVertexerComponent.cxx
  18     @author Sergey Gorbunov
  19     @brief  Component for reconstruct primary vertex and V0's
  20 */
  21
  22 #include "AliHLTGlobalVertexerComponent.h"
  23 #include "AliHLTDataTypes.h"
  24 #include "AliHLTExternalTrackParam.h"
  25 #include "AliHLTGlobalBarrelTrack.h"
  26 #include "AliCDBEntry.h"
  27 #include "AliCDBManager.h"
  28 #include <TFile.h>
  29 #include <TString.h>
  30 #include "TObjString.h"
  31 #include "TObjArray.h"
  32 #include "AliESDEvent.h"
  33 #include "AliESDtrack.h"
  34 #include "AliESDVertex.h"
  35 #include "AliESDv0.h"
  36 #include "AliHLTMessage.h"
  37 #include "TMath.h"
  38 #include "AliKFParticle.h"
  39 #include "AliKFVertex.h"
  40 #include "TStopwatch.h"
  41
  42 using namespace std;
  43
  44 /** ROOT macro for the implementation of ROOT specific class methods */
  45 ClassImp(AliHLTGlobalVertexerComponent)
  46
  47 AliHLTGlobalVertexerComponent::AliHLTGlobalVertexerComponent()
  48 :
  49   fNTracks(0),
  50   fTrackInfos(0),
  51   fPrimaryVtx(),
  52   fFitTracksToVertex(1),
  53   fConstrainedTrackDeviation(4.),
  54   fV0DaughterPrimDeviation( 2.5 ),
  55   fV0PrimDeviation( 3.5 ),
  56   fV0Chi(3.5),
  57   fV0DecayLengthInSigmas(3.),
  58   fV0TimeLimit(10.e-3),
  59   fBenchmark("GlobalVertexer")
  60 {
  61   // see header file for class documentation
  62   // or
  63   // refer to README to build package
  64   // or
  65   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  66
  67 }
  68
  69 AliHLTGlobalVertexerComponent::~AliHLTGlobalVertexerComponent()
  70 {
  71   // see header file for class documentation
  72
  73   if( fTrackInfos ) delete[] fTrackInfos;
  74 }
  75
  76 // Public functions to implement AliHLTComponent's interface.
  77 // These functions are required for the registration process
  78
  79 const char* AliHLTGlobalVertexerComponent::GetComponentID()
  80 {
  81   // see header file for class documentation
  82
  83   return "GlobalVertexer";
  84 }
  85
  86 void AliHLTGlobalVertexerComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
  87 {
  88   // see header file for class documentation
  89   list.clear();
  90   list.push_back( kAliHLTDataTypeESDObject );
  91   list.push_back( kAliHLTDataTypeTrack|kAliHLTDataOriginITS );
  92   list.push_back( kAliHLTDataTypeTrack|kAliHLTDataOriginTPC );
  93 }
  94
  95 AliHLTComponentDataType AliHLTGlobalVertexerComponent::GetOutputDataType()
  96 {
  97   // see header file for class documentation
  98   return kAliHLTMultipleDataType;
  99 }
 100
 101 int AliHLTGlobalVertexerComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
 102
 103 {
 104   // see header file for class documentation
 105   tgtList.clear();
 106   tgtList.push_back( kAliHLTDataTypeGlobalVertexer|kAliHLTDataOriginOut);
 107   tgtList.push_back( kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut);
 108   tgtList.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut);
 109   return tgtList.size();
 110 }
 111
 112 void AliHLTGlobalVertexerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 113 {
 114   // see header file for class documentation
 115   // XXX TODO: Find more realistic values.
 116   constBase = 80000;
 117   inputMultiplier = 2.;
 118 }
 119
 120 AliHLTComponent* AliHLTGlobalVertexerComponent::Spawn()
 121 {
 122   // see header file for class documentation
 123   return new AliHLTGlobalVertexerComponent;
 124 }
 125
 126 int AliHLTGlobalVertexerComponent::DoInit( int argc, const char** argv )
 127 {
 128   // init
 129
 130   fBenchmark.Reset();
 131   fBenchmark.SetTimer(0,"total");
 132   fBenchmark.SetTimer(1,"convert");
 133   fBenchmark.SetTimer(2,"vprim");
 134   fBenchmark.SetTimer(3,"v0");
 135   fV0TimeLimit = 10.e-3;
 136
 137   int iResult=0;
 138   TString configuration="";
 139   TString argument="";
 140   for (int i=0; i<argc && iResult>=0; i++) {
 141     argument=argv[i];
 142     if (!configuration.IsNull()) configuration+=" ";
 143     configuration+=argument;
 144   }
 145
 146   if (!configuration.IsNull()) {
 147     iResult=Configure(configuration.Data());
 148   }
 149
 150   return iResult;
 151 }
 152
 153 int AliHLTGlobalVertexerComponent::DoDeinit()
 154 {
 155   // see header file for class documentation
 156
 157   if( fTrackInfos ) delete[] fTrackInfos;
 158
 159   fTrackInfos = 0;
 160   fFitTracksToVertex = 1;
 161   fConstrainedTrackDeviation = 4.;
 162   fV0DaughterPrimDeviation = 2.5 ;
 163   fV0PrimDeviation =3.5;
 164   fV0Chi = 3.5;
 165   fV0DecayLengthInSigmas = 3.;
 166   fV0TimeLimit = 10.e-3;
 167
 168   return 0;
 169 }
 170
 171 int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*evtData*/,
 172                                             AliHLTComponentTriggerData& /*trigData*/ )
 173 {
 174   //cout<<"AliHLTGlobalVertexerComponent::DoEvent called"<<endl;
 175
 176   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
 177     return 0;
 178
 179   int iResult = 0;
 180   //cout<<"SG: GlobalVertexer:: DoEvent called"<<endl;
 181   fBenchmark.StartNewEvent();
 182   fBenchmark.Start(0);
 183
 184   vector< AliExternalTrackParam > tracks;
 185   vector< int > trackId;
 186   vector< pair<int,int> > v0s;
 187
 188   AliESDEvent *event = 0;
 189
 190   for( const AliHLTComponentBlockData *i= GetFirstInputBlock(kAliHLTDataTypeESDObject); i!=NULL; i=GetNextInputBlock() ){
 191     fBenchmark.AddInput(i->fSize);
 192   }
 193
 194
 195   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
 196     event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
 197     if( !event ) continue;
 198
 199     event->GetStdContent();
 200     Int_t nESDTracks=event->GetNumberOfTracks();
 201
 202     for (Int_t iTr=0; iTr<nESDTracks; iTr++){
 203       AliESDtrack *pTrack = event->GetTrack(iTr);
 204       if( !pTrack  ) continue;
 205       if( pTrack->GetKinkIndex(0)>0) continue;
 206       if( !( pTrack->GetStatus()&AliESDtrack::kTPCin ) ) continue;
 207       tracks.push_back(*pTrack);
 208       trackId.push_back(iTr);
 209     }
 210     break;
 211   }
 212
 213   if( tracks.size()==0 ){
 214
 215     for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
 216          pBlock!=NULL; pBlock=GetNextInputBlock()) {
 217
 218       fBenchmark.AddInput(pBlock->fSize);
 219
 220       AliHLTTracksData* dataPtr = reinterpret_cast<AliHLTTracksData*>( pBlock->fPtr );
 221       int nTracks = dataPtr->fCount;
 222       AliHLTExternalTrackParam* currOutTrack = dataPtr->fTracklets;
 223       for( int itr=0; itr<nTracks; itr++ ){
 224         AliHLTGlobalBarrelTrack t(*currOutTrack);
 225         tracks.push_back( t );
 226         //!SW trackId.push_back( currOutTrack->fTrackID );
 227         trackId.push_back( itr );
 228         unsigned int dSize = sizeof( AliHLTExternalTrackParam ) + currOutTrack->fNPoints * sizeof( unsigned int );
 229         currOutTrack = ( AliHLTExternalTrackParam* )( (( Byte_t * )currOutTrack) + dSize );
 230       }
 231     }
 232   }
 233
 234   if( tracks.size()==0 ){
 235     for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
 236          pBlock!=NULL; pBlock=GetNextInputBlock()) {
 237
 238       fBenchmark.AddInput(pBlock->fSize);
 239
 240       AliHLTTracksData* dataPtr = reinterpret_cast<AliHLTTracksData*>( pBlock->fPtr );
 241       int nTracks = dataPtr->fCount;
 242       AliHLTExternalTrackParam* currOutTrack = dataPtr->fTracklets;
 243       for( int itr=0; itr<nTracks; itr++ ){
 244         AliHLTGlobalBarrelTrack t(*currOutTrack);
 245         tracks.push_back( t );
 246         //!SW trackId.push_back( currOutTrack->fTrackID );
 247         trackId.push_back( itr );
 248         unsigned int dSize = sizeof( AliHLTExternalTrackParam ) + currOutTrack->fNPoints * sizeof( unsigned int );
 249         currOutTrack = ( AliHLTExternalTrackParam* )( (( Byte_t * )currOutTrack) + dSize );
 250       }
 251     }
 252   }
 253
 254
 255   // primary vertex & V0's
 256
 257   AliKFParticle::SetField( GetBz() );
 258
 259   fBenchmark.Start(1);
 260
 261   {  //* Fill fTrackInfo array
 262
 263     if( fTrackInfos ) delete[] fTrackInfos;
 264     fTrackInfos = 0;
 265     fNTracks=tracks.size();
 266     fTrackInfos = new AliESDTrackInfo[ fNTracks ];
 267     for (Int_t iTr=0; iTr<fNTracks; iTr++){
 268       AliESDTrackInfo &info = fTrackInfos[iTr];
 269       info.fOK = 1;
 270       info.fPrimUsedFlag = 0;
 271       info.fParticle = AliKFParticle( tracks[iTr], 211 );
 272       for( int i=0; i<8; i++ ) if( !finite(info.fParticle.GetParameter(i)) ) info.fOK = 0;
 273       for( int i=0; i<36; i++ ) if( !finite(info.fParticle.GetCovariance(i)) ) info.fOK = 0;
 274     }
 275   }
 276
 277   fBenchmark.Stop(1);
 278   fBenchmark.Start(2);
 279   FindPrimaryVertex();
 280   fBenchmark.Stop(2);
 281   fBenchmark.Start(3);
 282   FindV0s( v0s );
 283   fBenchmark.Stop(3);
 284
 285   int *buf = new int[sizeof(AliHLTGlobalVertexerData)/sizeof(int)+1 + fNTracks + 2*v0s.size()];
 286   AliHLTGlobalVertexerData *data = reinterpret_cast<AliHLTGlobalVertexerData*>(buf);
 287
 288   if( data) {  // fill the output structure
 289
 290     data->fFitTracksToVertex = fFitTracksToVertex;
 291     for( int i=0; i<3; i++ ) data->fPrimP[i] = fPrimaryVtx.Parameters()[i];
 292     for( int i=0; i<6; i++ ) data->fPrimC[i] = fPrimaryVtx.CovarianceMatrix()[i];
 293     data->fPrimChi2 = fPrimaryVtx.GetChi2();
 294     data->fPrimNContributors = fPrimaryVtx.GetNContributors();
 295     data->fNPrimTracks = 0;
 296     for( Int_t i = 0; i<fNTracks; i++ ){
 297       if( !fTrackInfos[i].fPrimUsedFlag ) continue;
 298       if( fTrackInfos[i].fPrimDeviation > fConstrainedTrackDeviation ) continue;
 299       data->fTrackIndices[ (data->fNPrimTracks)++ ] =  trackId[i];
 300     }
 301     int *listV0 = data->fTrackIndices + data->fNPrimTracks;
 302     data->fNV0s = v0s.size();
 303     for( int i=0; i<data->fNV0s; i++ ){
 304       listV0[2*i] = trackId[v0s[i].first];
 305       listV0[2*i+1] = trackId[v0s[i].second];
 306     }
 307
 308     unsigned int blockSize = sizeof(AliHLTGlobalVertexerData) + (data->fNPrimTracks + 2*data->fNV0s)*sizeof(int);
 309
 310     iResult = PushBack( reinterpret_cast<void*>(data), blockSize, kAliHLTDataTypeGlobalVertexer|kAliHLTDataOriginOut );
 311     fBenchmark.AddOutput(blockSize);
 312   }
 313
 314
 315   // output the vertex if found
 316   {
 317     if( iResult==0 && data && data->fPrimNContributors >=3 ){
 318       AliESDVertex vESD( data->fPrimP, data->fPrimC, data->fPrimChi2, data->fPrimNContributors );
 319       iResult = PushBack( (TObject*) &vESD, kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut,0 );
 320       fBenchmark.AddOutput(GetLastObjectSize());
 321     }
 322   }
 323
 324   // output the ESD event
 325   if( iResult==0 && event && data ){
 326     FillESD( event, data );
 327     iResult = PushBack( event, kAliHLTDataTypeESDObject|kAliHLTDataOriginOut, 0);
 328     fBenchmark.AddOutput(GetLastObjectSize());
 329   }
 330
 331   delete[] buf;
 332
 333   fBenchmark.Stop(0);
 334   HLTInfo(fBenchmark.GetStatistics());
 335   return iResult;
 336 }
 337
 338
 339 int AliHLTGlobalVertexerComponent::Configure(const char* arguments)
 340 {
 341
 342   int iResult=0;
 343   if (!arguments) return iResult;
 344
 345   TString allArgs=arguments;
 346   TString argument;
 347
 348   TObjArray* pTokens=allArgs.Tokenize(" ");
 349   int bMissingParam=0;
 350
 351   if (pTokens) {
 352     for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
 353       argument=((TObjString*)pTokens->At(i))->GetString();
 354       if (argument.IsNull()) continue;
 355
 356       if (argument.CompareTo("-fitTracksToVertex")==0) {
 357         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 358         HLTInfo("fitTracksToVertex is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 359         fFitTracksToVertex=((TObjString*)pTokens->At(i))->GetString().Atoi();
 360         continue;
 361       }
 362       else if (argument.CompareTo("-constrainedTrackDeviation")==0) {
 363         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 364         HLTInfo("constrainedTrackDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 365         fConstrainedTrackDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 366         continue;
 367       }
 368       else if (argument.CompareTo("-v0DaughterPrimDeviation")==0) {
 369         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 370         HLTInfo("v0DaughterPrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 371         fV0DaughterPrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 372         continue;
 373       }
 374       else if (argument.CompareTo("-v0PrimDeviation")==0) {
 375         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 376         HLTInfo("v0PrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 377         fV0PrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 378         continue;
 379       }
 380       else if (argument.CompareTo("-v0Chi")==0) {
 381         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 382         HLTInfo("v0Chi is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 383         fV0Chi=((TObjString*)pTokens->At(i))->GetString().Atof();
 384         continue;
 385       }
 386       else if (argument.CompareTo("-v0DecayLengthInSigmas")==0) {
 387         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 388         HLTInfo("v0DecayLengthInSigmas is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 389         fV0DecayLengthInSigmas=((TObjString*)pTokens->At(i))->GetString().Atof();
 390         continue;
 391       }
 392       else if (argument.CompareTo("-v0MinEventRate")==0) {
 393         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 394         HLTInfo("Minimum event rate for V0 finder is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 395         Double_t rate = ((TObjString*)pTokens->At(i))->GetString().Atof();
 396         fV0TimeLimit = (rate >0 ) ?1./rate :60; // 1 minute maximum time
 397         continue;
 398       }
 399       else {
 400         HLTError("unknown argument %s", argument.Data());
 401         iResult=-EINVAL;
 402         break;
 403       }
 404     }
 405     delete pTokens;
 406   }
 407   if (bMissingParam) {
 408     HLTError("missing parameter for argument %s", argument.Data());
 409     iResult=-EINVAL;
 410   }
 411
 412   return iResult;
 413 }
 414
 415 int AliHLTGlobalVertexerComponent::Reconfigure(const char* /*cdbEntry*/, const char* /*chainId*/)
 416 {
 417   // see header file for class documentation
 418
 419   return 0; // no CDB path is set so far
 420   /*
 421   int iResult=0;
 422   const char* path="HLT/ConfigTPC/KryptonHistoComponent";
 423   const char* defaultNotify="";
 424   if (cdbEntry) {
 425     path=cdbEntry;
 426     defaultNotify=" (default)";
 427   }
 428   if (path) {
 429     HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
 430     AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path);
 431     if (pEntry) {
 432       TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
 433       if (pString) {
 434         HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
 435         iResult=Configure(pString->GetString().Data());
 436       } else {
 437         HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
 438       }
 439     } else {
 440       HLTError("can not fetch object \"%s\" from CDB", path);
 441     }
 442   }
 443
 444   return iResult;
 445 */
 446 }
 447
 448
 449 struct AliHLTGlobalVertexerDeviation
 450 {
 451   int fI; // track index
 452   float fD; // deviation from vertex
 453
 454   bool operator<(const AliHLTGlobalVertexerDeviation &a) const { return fD<a.fD; }
 455 };
 456
 457
 458 void AliHLTGlobalVertexerComponent::FindPrimaryVertex()
 459 {
 460   //* Find event primary vertex
 461
 462   fPrimaryVtx.Initialize();
 463   //fPrimaryVtx.SetBeamConstraint(fESD->GetDiamondX(),fESD->GetDiamondY(),0,
 464   //TMath::Sqrt(fESD->GetSigma2DiamondX()),TMath::Sqrt(fESD->GetSigma2DiamondY()),5.3);
 465
 466   // select rough region (in sigmas) in which the vertex could be found, all tracks outside these limits are rejected
 467   // from the primary vertex finding
 468   fPrimaryVtx.SetBeamConstraint( 0, 0, 0, 3., 3., 5.3 );
 469
 470   const AliKFParticle **vSelected = new const AliKFParticle*[fNTracks]; //* Selected particles for the vertex fit
 471   AliHLTGlobalVertexerDeviation *dev = new AliHLTGlobalVertexerDeviation[fNTracks];
 472
 473   Int_t nSelected = 0;
 474
 475   for( Int_t i = 0; i<fNTracks; i++){
 476     if(!fTrackInfos[i].fOK ) continue;
 477     //if( fESD->GetTrack(i)->GetTPCNcls()<60  ) continue;
 478     const AliKFParticle &p = fTrackInfos[i].fParticle;
 479     Double_t chi = p.GetDeviationFromVertex( fPrimaryVtx );
 480     if( chi > fConstrainedTrackDeviation ) continue;
 481     dev[nSelected].fI = i;
 482     dev[nSelected].fD = chi;
 483     vSelected[nSelected] = &(fTrackInfos[i].fParticle);
 484     nSelected++;
 485   }
 486
 487   // fit
 488
 489   while( nSelected>2 ){
 490
 491     //* Primary vertex finder with rejection of outliers
 492
 493     for( Int_t i = 0; i<nSelected; i++){
 494       vSelected[i] = &(fTrackInfos[dev[i].fI].fParticle);
 495     }
 496
 497     double xv = fPrimaryVtx.GetX();
 498     double yv = fPrimaryVtx.GetY();
 499     double zv = fPrimaryVtx.GetZ(); // values from previous iteration of calculations
 500     fPrimaryVtx.Initialize();
 501     fPrimaryVtx.SetBeamConstraint( 0, 0, 0, 3., 3., 5.3 );
 502     fPrimaryVtx.SetVtxGuess( xv, yv, zv );
 503
 504     fPrimaryVtx.Construct( vSelected, nSelected, 0, -1, 1 ); // refilled for every iteration
 505
 506     for( Int_t it=0; it<nSelected; it++ ){
 507       const AliKFParticle &p = fTrackInfos[dev[it].fI].fParticle;
 508       if( nSelected <= 20 ){
 509         AliKFVertex tmp =  fPrimaryVtx - p; // exclude the current track from the sample and recalculate the vertex
 510         dev[it].fD = p.GetDeviationFromVertex( tmp );
 511       } else {
 512         dev[it].fD = p.GetDeviationFromVertex( fPrimaryVtx );
 513       }
 514     }
 515     sort(dev,dev+nSelected); // sort tracks with increasing chi2 (used for rejection)
 516
 517     int nRemove = (int) ( 0.3*nSelected );  //remove 30% of the tracks (done for performance, only if there are more than 20 tracks)
 518     if( nSelected - nRemove <=20 ) nRemove = 1;  // removal based on the chi2 of every track
 519     int firstRemove = nSelected - nRemove;
 520     while( firstRemove<nSelected ){
 521       if( dev[firstRemove].fD >= fConstrainedTrackDeviation ) break;
 522       firstRemove++;
 523     }
 524     if( firstRemove>=nSelected ) break;
 525     nSelected = firstRemove;
 526   }
 527
 528   for( Int_t i = 0; i<fNTracks; i++){
 529     fTrackInfos[i].fPrimUsedFlag = 0;
 530   }
 531
 532   if( nSelected < 3 ){  // no vertex for fewer than 3 contributors
 533     fPrimaryVtx.NDF() = -3;
 534     fPrimaryVtx.Chi2() = 0;
 535     nSelected = 0;
 536   }
 537
 538   for( Int_t i = 0; i<nSelected; i++){
 539     AliESDTrackInfo &info = fTrackInfos[dev[i].fI];
 540     info.fPrimUsedFlag = 1;
 541     info.fPrimDeviation = dev[i].fD;
 542   }
 543
 544   for( Int_t i = 0; i<fNTracks; i++ ){
 545     AliESDTrackInfo &info = fTrackInfos[i];
 546     if( info.fPrimUsedFlag ) continue;
 547     info.fPrimDeviation = info.fParticle.GetDeviationFromVertex( fPrimaryVtx );
 548   }
 549
 550   delete[] vSelected;
 551   delete[] dev;
 552 }
 553
 554
 555
 556 void AliHLTGlobalVertexerComponent::FindV0s( vector<pair<int,int> > &v0s  )
 557 {
 558   //* V0 finder
 559
 560   AliKFVertex &primVtx = fPrimaryVtx;
 561   if( primVtx.GetNContributors()<3 ) return;
 562
 563   TStopwatch timer;
 564   Int_t statN = 0;
 565   Bool_t run = 1;
 566
 567   for( Int_t iTr = 0; iTr<fNTracks && run; iTr++ ){ //* first daughter
 568
 569     AliESDTrackInfo &info = fTrackInfos[iTr];
 570     if( !info.fOK ) continue;
 571     if( info.fParticle.GetQ() >0 ) continue;
 572     if( info.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
 573
 574     for( Int_t jTr = 0; jTr<fNTracks; jTr++ ){  //* second daughter
 575
 576
 577       AliESDTrackInfo &jnfo = fTrackInfos[jTr];
 578       if( !jnfo.fOK ) continue;
 579       if( jnfo.fParticle.GetQ() < 0 ) continue;
 580       if( jnfo.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
 581
 582       // check the time once a while...
 583
 584       if( (++statN)%100 ==0 ){
 585         if( timer.RealTime()>= fV0TimeLimit ){  run = 0; break; }
 586         timer.Continue();
 587       }
 588
 589       //* check if the particles fit
 590
 591       if( info.fParticle.GetDeviationFromParticle(jnfo.fParticle) > fV0Chi ) continue;
 592
 593       //* construct V0 mother
 594
 595       AliKFParticle v0( info.fParticle, jnfo.fParticle );
 596
 597       //* check V0 Chi^2
 598
 599       if( v0.GetChi2()<0 || v0.GetChi2() > fV0Chi*fV0Chi*v0.GetNDF() ) continue;
 600
 601       //* subtruct daughters from primary vertex
 602
 603       AliKFVertex primVtxCopy = primVtx;
 604
 605       if( info.fPrimUsedFlag ){
 606         if( primVtxCopy.GetNContributors()<=2 ) continue;
 607         primVtxCopy -= info.fParticle;
 608       }
 609       if( jnfo.fPrimUsedFlag ){
 610         if( primVtxCopy.GetNContributors()<=2 ) continue;
 611         primVtxCopy -= jnfo.fParticle;
 612       }
 613       //* Check v0 Chi^2 deviation from primary vertex
 614
 615       if( v0.GetDeviationFromVertex( primVtxCopy ) > fV0PrimDeviation ) continue;
 616
 617       //* Add V0 to primary vertex to improve the primary vertex resolution
 618
 619       primVtxCopy += v0;
 620
 621       //* Set production vertex for V0
 622
 623       v0.SetProductionVertex( primVtxCopy );
 624
 625       //* Get V0 decay length with estimated error
 626
 627       Double_t length, sigmaLength;
 628       if( v0.GetDecayLength( length, sigmaLength ) ) continue;
 629
 630       //* Reject V0 if it decays too close[sigma] to the primary vertex
 631
 632       if( length  < fV0DecayLengthInSigmas*sigmaLength ) continue;
 633
 634       //* keep v0
 635
 636       v0s.push_back(pair<int,int>(iTr,jTr));
 637     }
 638   }
 639 }
 640
 641
 642
 643
 644 void AliHLTGlobalVertexerComponent::FillESD( AliESDEvent *event, AliHLTGlobalVertexerData *data
 645 )
 646 {
 647   //* put output of a vertexer to the esd event
 648
 649   Int_t nESDTracks = event->GetNumberOfTracks();
 650
 651   const int *listPrim = data->fTrackIndices;
 652   const int *listV0 = data->fTrackIndices + data->fNPrimTracks;
 653
 654   std::map<int,int> mapId;
 655   bool *constrainedToVtx   = new bool[nESDTracks];
 656
 657   for( int i=0; i<nESDTracks; i++ ){
 658     constrainedToVtx[i] = 0;
 659     if( !event->GetTrack(i) ) continue;
 660     mapId[ event->GetTrack(i)->GetID() ] = i;
 661   }
 662
 663   if( data->fPrimNContributors >=3 ){
 664
 665     AliESDVertex vESD( data->fPrimP, data->fPrimC, data->fPrimChi2, data->fPrimNContributors );
 666     event->SetPrimaryVertexTPC( &vESD );
 667     event->SetPrimaryVertexTracks( &vESD );
 668
 669     // relate tracks to the primary vertex
 670
 671     if( data->fFitTracksToVertex ){
 672       for( Int_t i = 0; i<data->fNPrimTracks; i++ ){
 673         Int_t id = listPrim[ i ];
 674         map<int,int>::iterator it = mapId.find(id);
 675         if( it==mapId.end() ) continue;
 676         Int_t itr = it->second;
 677         event->GetTrack(itr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
 678         constrainedToVtx[ itr ] = 1;
 679       }
 680     }
 681   }
 682
 683   //* add ESD v0s and relate tracks to v0s
 684
 685
 686   for( int i=0; i<data->fNV0s; i++ ){
 687
 688     Int_t id1 = listV0[ 2*i ];
 689     Int_t id2 = listV0[ 2*i + 1];
 690     map<int,int>::iterator it = mapId.find(id1);
 691     if( it==mapId.end() ) continue;
 692     Int_t iTr = it->second;
 693     it = mapId.find(id2);
 694     if( it==mapId.end() ) continue;
 695     Int_t jTr = it->second;
 696
 697     AliESDv0 v0( *event->GetTrack( iTr ), iTr, *event->GetTrack( jTr ), jTr );
 698     event->AddV0( &v0 );
 699
 700     // relate the tracks to the vertex
 701
 702     if( data->fFitTracksToVertex ){
 703       if( constrainedToVtx[iTr] || constrainedToVtx[jTr] ) continue;
 704       double pos[3];
 705       double sigma[3] = {.1,.1,.1};
 706       v0.XvYvZv(pos);
 707       AliESDVertex vESD(pos, sigma);
 708       event->GetTrack(iTr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
 709       event->GetTrack(jTr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
 710       constrainedToVtx[iTr] = 1;
 711       constrainedToVtx[jTr] = 1;
 712     }
 713   }
 714
 715   delete[] constrainedToVtx;
 716 }