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
 223       AliHLTExternalTrackParam* currOutTrack = dataPtr->fTracklets;
 224       for( int itr=0; itr<nTracks; itr++ ){
 225         AliHLTGlobalBarrelTrack t(*currOutTrack);
 226         tracks.push_back( t );
 227         trackId.push_back( currOutTrack->fTrackID );
 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         trackId.push_back( currOutTrack->fTrackID );
 247         unsigned int dSize = sizeof( AliHLTExternalTrackParam ) + currOutTrack->fNPoints * sizeof( unsigned int );
 248         currOutTrack = ( AliHLTExternalTrackParam* )( (( Byte_t * )currOutTrack) + dSize );
 249       }
 250     }
 251   }
 252
 253
 254   // primary vertex & V0's
 255
 256   AliKFParticle::SetField( GetBz() );
 257
 258   fBenchmark.Start(1);
 259
 260   {  //* Fill fTrackInfo array
 261
 262     if( fTrackInfos ) delete[] fTrackInfos;
 263     fTrackInfos = 0;
 264     fNTracks=tracks.size();
 265     fTrackInfos = new AliESDTrackInfo[ fNTracks ];
 266     for (Int_t iTr=0; iTr<fNTracks; iTr++){
 267       AliESDTrackInfo &info = fTrackInfos[iTr];
 268       info.fOK = 1;
 269       info.fPrimUsedFlag = 0;
 270       info.fParticle = AliKFParticle( tracks[iTr], 211 );
 271       for( int i=0; i<8; i++ ) if( !finite(info.fParticle.GetParameter(i)) ) info.fOK = 0;
 272       for( int i=0; i<36; i++ ) if( !finite(info.fParticle.GetCovariance(i)) ) info.fOK = 0;
 273     }
 274   }
 275
 276   fBenchmark.Stop(1);
 277   fBenchmark.Start(2);
 278   FindPrimaryVertex();
 279   fBenchmark.Stop(2);
 280   fBenchmark.Start(3);
 281   FindV0s( v0s );
 282   fBenchmark.Stop(3);
 283
 284   int *buf = new int[sizeof(AliHLTGlobalVertexerData)/sizeof(int)+1 + fNTracks + 2*v0s.size()];
 285   AliHLTGlobalVertexerData *data = reinterpret_cast<AliHLTGlobalVertexerData*>(buf);
 286
 287   if( data) {  // fill the output structure
 288
 289     data->fFitTracksToVertex = fFitTracksToVertex;
 290     for( int i=0; i<3; i++ ) data->fPrimP[i] = fPrimaryVtx.Parameters()[i];
 291     for( int i=0; i<6; i++ ) data->fPrimC[i] = fPrimaryVtx.CovarianceMatrix()[i];
 292     data->fPrimChi2 = fPrimaryVtx.GetChi2();
 293     data->fPrimNContributors = fPrimaryVtx.GetNContributors();
 294     data->fNPrimTracks = 0;
 295     for( Int_t i = 0; i<fNTracks; i++ ){
 296       if( !fTrackInfos[i].fPrimUsedFlag ) continue;
 297       if( fTrackInfos[i].fPrimDeviation > fConstrainedTrackDeviation ) continue;
 298       data->fTrackIndices[ (data->fNPrimTracks)++ ] =  trackId[i];
 299     }
 300     int *listV0 = data->fTrackIndices + data->fNPrimTracks;
 301     data->fNV0s = v0s.size();
 302     for( int i=0; i<data->fNV0s; i++ ){
 303       listV0[2*i] = trackId[v0s[i].first];
 304       listV0[2*i+1] = trackId[v0s[i].second];
 305     }
 306
 307     unsigned int blockSize = sizeof(AliHLTGlobalVertexerData) + (data->fNPrimTracks + 2*data->fNV0s)*sizeof(int);
 308
 309     iResult = PushBack( reinterpret_cast<void*>(data), blockSize, kAliHLTDataTypeGlobalVertexer|kAliHLTDataOriginOut );
 310     fBenchmark.AddOutput(blockSize);
 311   }
 312
 313
 314   // output the vertex if found
 315   {
 316     if( iResult==0 && data && data->fPrimNContributors >=3 ){
 317       AliESDVertex vESD( data->fPrimP, data->fPrimC, data->fPrimChi2, data->fPrimNContributors );
 318       iResult = PushBack( (TObject*) &vESD, kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut,0 );
 319       fBenchmark.AddOutput(GetLastObjectSize());
 320     }
 321   }
 322
 323   // output the ESD event
 324   if( iResult==0 && event && data ){
 325     FillESD( event, data );
 326     iResult = PushBack( event, kAliHLTDataTypeESDObject|kAliHLTDataOriginOut, 0);
 327     fBenchmark.AddOutput(GetLastObjectSize());
 328   }
 329
 330   delete[] buf;
 331
 332   fBenchmark.Stop(0);
 333   HLTInfo(fBenchmark.GetStatistics());
 334   return iResult;
 335 }
 336
 337
 338 int AliHLTGlobalVertexerComponent::Configure(const char* arguments)
 339 {
 340
 341   int iResult=0;
 342   if (!arguments) return iResult;
 343
 344   TString allArgs=arguments;
 345   TString argument;
 346
 347   TObjArray* pTokens=allArgs.Tokenize(" ");
 348   int bMissingParam=0;
 349
 350   if (pTokens) {
 351     for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
 352       argument=((TObjString*)pTokens->At(i))->GetString();
 353       if (argument.IsNull()) continue;
 354
 355       if (argument.CompareTo("-fitTracksToVertex")==0) {
 356         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 357         HLTInfo("fitTracksToVertex is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 358         fFitTracksToVertex=((TObjString*)pTokens->At(i))->GetString().Atoi();
 359         continue;
 360       }
 361       else if (argument.CompareTo("-constrainedTrackDeviation")==0) {
 362         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 363         HLTInfo("constrainedTrackDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 364         fConstrainedTrackDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 365         continue;
 366       }
 367       else if (argument.CompareTo("-v0DaughterPrimDeviation")==0) {
 368         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 369         HLTInfo("v0DaughterPrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 370         fV0DaughterPrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 371         continue;
 372       }
 373       else if (argument.CompareTo("-v0PrimDeviation")==0) {
 374         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 375         HLTInfo("v0PrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 376         fV0PrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 377         continue;
 378       }
 379       else if (argument.CompareTo("-v0Chi")==0) {
 380         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 381         HLTInfo("v0Chi is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 382         fV0Chi=((TObjString*)pTokens->At(i))->GetString().Atof();
 383         continue;
 384       }
 385       else if (argument.CompareTo("-v0DecayLengthInSigmas")==0) {
 386         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 387         HLTInfo("v0DecayLengthInSigmas is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 388         fV0DecayLengthInSigmas=((TObjString*)pTokens->At(i))->GetString().Atof();
 389         continue;
 390       }
 391       else if (argument.CompareTo("-v0MinEventRate")==0) {
 392         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 393         HLTInfo("Minimum event rate for V0 finder is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 394         Double_t rate = ((TObjString*)pTokens->At(i))->GetString().Atof();
 395         fV0TimeLimit = (rate >0 ) ?1./rate :60; // 1 minute maximum time
 396         continue;
 397       }
 398       else {
 399         HLTError("unknown argument %s", argument.Data());
 400         iResult=-EINVAL;
 401         break;
 402       }
 403     }
 404     delete pTokens;
 405   }
 406   if (bMissingParam) {
 407     HLTError("missing parameter for argument %s", argument.Data());
 408     iResult=-EINVAL;
 409   }
 410
 411   return iResult;
 412 }
 413
 414 int AliHLTGlobalVertexerComponent::Reconfigure(const char* /*cdbEntry*/, const char* /*chainId*/)
 415 {
 416   // see header file for class documentation
 417
 418   return 0; // no CDB path is set so far
 419   /*
 420   int iResult=0;
 421   const char* path="HLT/ConfigTPC/KryptonHistoComponent";
 422   const char* defaultNotify="";
 423   if (cdbEntry) {
 424     path=cdbEntry;
 425     defaultNotify=" (default)";
 426   }
 427   if (path) {
 428     HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
 429     AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path);
 430     if (pEntry) {
 431       TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
 432       if (pString) {
 433         HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
 434         iResult=Configure(pString->GetString().Data());
 435       } else {
 436         HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
 437       }
 438     } else {
 439       HLTError("can not fetch object \"%s\" from CDB", path);
 440     }
 441   }
 442
 443   return iResult;
 444 */
 445 }
 446
 447
 448 struct AliHLTGlobalVertexerDeviation
 449 {
 450   int fI; // track index
 451   float fD; // deviation from vertex
 452
 453   bool operator<(const AliHLTGlobalVertexerDeviation &a) const { return fD<a.fD; }
 454 };
 455
 456
 457 void AliHLTGlobalVertexerComponent::FindPrimaryVertex()
 458 {
 459   //* Find event primary vertex
 460
 461   fPrimaryVtx.Initialize();
 462   //fPrimaryVtx.SetBeamConstraint(fESD->GetDiamondX(),fESD->GetDiamondY(),0,
 463   //TMath::Sqrt(fESD->GetSigma2DiamondX()),TMath::Sqrt(fESD->GetSigma2DiamondY()),5.3);
 464
 465   // select rough region (in sigmas) in which the vertex could be found, all tracks outside these limits are rejected
 466   // from the primary vertex finding
 467   fPrimaryVtx.SetBeamConstraint( 0, 0, 0, 3., 3., 5.3 );
 468
 469   const AliKFParticle **vSelected = new const AliKFParticle*[fNTracks]; //* Selected particles for the vertex fit
 470   AliHLTGlobalVertexerDeviation *dev = new AliHLTGlobalVertexerDeviation[fNTracks];
 471
 472   Int_t nSelected = 0;
 473
 474   for( Int_t i = 0; i<fNTracks; i++){
 475     if(!fTrackInfos[i].fOK ) continue;
 476     //if( fESD->GetTrack(i)->GetTPCNcls()<60  ) continue;
 477     const AliKFParticle &p = fTrackInfos[i].fParticle;
 478     Double_t chi = p.GetDeviationFromVertex( fPrimaryVtx );
 479     if( chi > fConstrainedTrackDeviation ) continue;
 480     dev[nSelected].fI = i;
 481     dev[nSelected].fD = chi;
 482     vSelected[nSelected] = &(fTrackInfos[i].fParticle);
 483     nSelected++;
 484   }
 485
 486   // fit
 487
 488   while( nSelected>2 ){
 489
 490     //* Primary vertex finder with rejection of outliers
 491
 492     for( Int_t i = 0; i<nSelected; i++){
 493       vSelected[i] = &(fTrackInfos[dev[i].fI].fParticle);
 494     }
 495
 496     double xv = fPrimaryVtx.GetX();
 497     double yv = fPrimaryVtx.GetY();
 498     double zv = fPrimaryVtx.GetZ(); // values from previous iteration of calculations
 499     fPrimaryVtx.Initialize();
 500     fPrimaryVtx.SetBeamConstraint( 0, 0, 0, 3., 3., 5.3 );
 501     fPrimaryVtx.SetVtxGuess( xv, yv, zv );
 502
 503     fPrimaryVtx.Construct( vSelected, nSelected, 0, -1, 1 ); // refilled for every iteration
 504
 505     for( Int_t it=0; it<nSelected; it++ ){
 506       const AliKFParticle &p = fTrackInfos[dev[it].fI].fParticle;
 507       if( nSelected <= 20 ){
 508         AliKFVertex tmp =  fPrimaryVtx - p; // exclude the current track from the sample and recalculate the vertex
 509         dev[it].fD = p.GetDeviationFromVertex( tmp );
 510       } else {
 511         dev[it].fD = p.GetDeviationFromVertex( fPrimaryVtx );
 512       }
 513     }
 514     sort(dev,dev+nSelected); // sort tracks with increasing chi2 (used for rejection)
 515
 516     int nRemove = (int) ( 0.3*nSelected );  //remove 30% of the tracks (done for performance, only if there are more than 20 tracks)
 517     if( nSelected - nRemove <=20 ) nRemove = 1;  // removal based on the chi2 of every track
 518     int firstRemove = nSelected - nRemove;
 519     while( firstRemove<nSelected ){
 520       if( dev[firstRemove].fD >= fConstrainedTrackDeviation ) break;
 521       firstRemove++;
 522     }
 523     if( firstRemove>=nSelected ) break;
 524     nSelected = firstRemove;
 525   }
 526
 527   for( Int_t i = 0; i<fNTracks; i++){
 528     fTrackInfos[i].fPrimUsedFlag = 0;
 529   }
 530
 531   if( nSelected < 3 ){  // no vertex for fewer than 3 contributors
 532     fPrimaryVtx.NDF() = -3;
 533     fPrimaryVtx.Chi2() = 0;
 534     nSelected = 0;
 535   }
 536
 537   for( Int_t i = 0; i<nSelected; i++){
 538     AliESDTrackInfo &info = fTrackInfos[dev[i].fI];
 539     info.fPrimUsedFlag = 1;
 540     info.fPrimDeviation = dev[i].fD;
 541   }
 542
 543   for( Int_t i = 0; i<fNTracks; i++ ){
 544     AliESDTrackInfo &info = fTrackInfos[i];
 545     if( info.fPrimUsedFlag ) continue;
 546     info.fPrimDeviation = info.fParticle.GetDeviationFromVertex( fPrimaryVtx );
 547   }
 548
 549   delete[] vSelected;
 550   delete[] dev;
 551 }
 552
 553
 554
 555 void AliHLTGlobalVertexerComponent::FindV0s( vector<pair<int,int> > &v0s  )
 556 {
 557   //* V0 finder
 558
 559   AliKFVertex &primVtx = fPrimaryVtx;
 560   if( primVtx.GetNContributors()<3 ) return;
 561
 562   TStopwatch timer;
 563   Int_t statN = 0;
 564   Bool_t run = 1;
 565
 566   for( Int_t iTr = 0; iTr<fNTracks && run; iTr++ ){ //* first daughter
 567
 568     AliESDTrackInfo &info = fTrackInfos[iTr];
 569     if( !info.fOK ) continue;
 570     if( info.fParticle.GetQ() >0 ) continue;
 571     if( info.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
 572
 573     for( Int_t jTr = 0; jTr<fNTracks; jTr++ ){  //* second daughter
 574
 575
 576       AliESDTrackInfo &jnfo = fTrackInfos[jTr];
 577       if( !jnfo.fOK ) continue;
 578       if( jnfo.fParticle.GetQ() < 0 ) continue;
 579       if( jnfo.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
 580
 581       // check the time once a while...
 582
 583       if( (++statN)%100 ==0 ){
 584         if( timer.RealTime()>= fV0TimeLimit ){  run = 0; break; }
 585         timer.Continue();
 586       }
 587
 588       //* check if the particles fit
 589
 590       if( info.fParticle.GetDeviationFromParticle(jnfo.fParticle) > fV0Chi ) continue;
 591
 592       //* construct V0 mother
 593
 594       AliKFParticle v0( info.fParticle, jnfo.fParticle );
 595
 596       //* check V0 Chi^2
 597
 598       if( v0.GetChi2()<0 || v0.GetChi2() > fV0Chi*fV0Chi*v0.GetNDF() ) continue;
 599
 600       //* subtruct daughters from primary vertex
 601
 602       AliKFVertex primVtxCopy = primVtx;
 603
 604       if( info.fPrimUsedFlag ){
 605         if( primVtxCopy.GetNContributors()<=2 ) continue;
 606         primVtxCopy -= info.fParticle;
 607       }
 608       if( jnfo.fPrimUsedFlag ){
 609         if( primVtxCopy.GetNContributors()<=2 ) continue;
 610         primVtxCopy -= jnfo.fParticle;
 611       }
 612       //* Check v0 Chi^2 deviation from primary vertex
 613
 614       if( v0.GetDeviationFromVertex( primVtxCopy ) > fV0PrimDeviation ) continue;
 615
 616       //* Add V0 to primary vertex to improve the primary vertex resolution
 617
 618       primVtxCopy += v0;
 619
 620       //* Set production vertex for V0
 621
 622       v0.SetProductionVertex( primVtxCopy );
 623
 624       //* Get V0 decay length with estimated error
 625
 626       Double_t length, sigmaLength;
 627       if( v0.GetDecayLength( length, sigmaLength ) ) continue;
 628
 629       //* Reject V0 if it decays too close[sigma] to the primary vertex
 630
 631       if( length  < fV0DecayLengthInSigmas*sigmaLength ) continue;
 632
 633       //* keep v0
 634
 635       v0s.push_back(pair<int,int>(iTr,jTr));
 636     }
 637   }
 638 }
 639
 640
 641
 642
 643 void AliHLTGlobalVertexerComponent::FillESD( AliESDEvent *event, AliHLTGlobalVertexerData *data
 644 )
 645 {
 646   //* put output of a vertexer to the esd event
 647
 648   Int_t nESDTracks = event->GetNumberOfTracks();
 649
 650   const int *listPrim = data->fTrackIndices;
 651   const int *listV0 = data->fTrackIndices + data->fNPrimTracks;
 652
 653   std::map<int,int> mapId;
 654   bool *constrainedToVtx   = new bool[nESDTracks];
 655
 656   for( int i=0; i<nESDTracks; i++ ){
 657     constrainedToVtx[i] = 0;
 658     if( !event->GetTrack(i) ) continue;
 659     mapId[ event->GetTrack(i)->GetID() ] = i;
 660   }
 661
 662   if( data->fPrimNContributors >=3 ){
 663
 664     AliESDVertex vESD( data->fPrimP, data->fPrimC, data->fPrimChi2, data->fPrimNContributors );
 665     event->SetPrimaryVertexTPC( &vESD );
 666     event->SetPrimaryVertexTracks( &vESD );
 667
 668     // relate tracks to the primary vertex
 669
 670     if( data->fFitTracksToVertex ){
 671       for( Int_t i = 0; i<data->fNPrimTracks; i++ ){
 672         Int_t id = listPrim[ i ];
 673         map<int,int>::iterator it = mapId.find(id);
 674         if( it==mapId.end() ) continue;
 675         Int_t itr = it->second;
 676         event->GetTrack(itr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
 677         constrainedToVtx[ itr ] = 1;
 678       }
 679     }
 680   }
 681
 682   //* add ESD v0s and relate tracks to v0s
 683
 684
 685   for( int i=0; i<data->fNV0s; i++ ){
 686
 687     Int_t id1 = listV0[ 2*i ];
 688     Int_t id2 = listV0[ 2*i + 1];
 689     map<int,int>::iterator it = mapId.find(id1);
 690     if( it==mapId.end() ) continue;
 691     Int_t iTr = it->second;
 692     it = mapId.find(id2);
 693     if( it==mapId.end() ) continue;
 694     Int_t jTr = it->second;
 695
 696     AliESDv0 v0( *event->GetTrack( iTr ), iTr, *event->GetTrack( jTr ), jTr );
 697     event->AddV0( &v0 );
 698
 699     // relate the tracks to the vertex
 700
 701     if( data->fFitTracksToVertex ){
 702       if( constrainedToVtx[iTr] || constrainedToVtx[jTr] ) continue;
 703       double pos[3];
 704       double sigma[3] = {.1,.1,.1};
 705       v0.XvYvZv(pos);
 706       AliESDVertex vESD(pos, sigma);
 707       event->GetTrack(iTr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
 708       event->GetTrack(jTr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
 709       constrainedToVtx[iTr] = 1;
 710       constrainedToVtx[jTr] = 1;
 711     }
 712   }
 713
 714   delete[] constrainedToVtx;
 715 }