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