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 "AliCDBEntry.h"
  28 #include "AliCDBManager.h"
  29 #include <TFile.h>
  30 #include <TString.h>
  31 #include "TObjString.h"
  32 #include "TObjArray.h"
  33 #include "TH1F.h"
  34 #include "TH2F.h"
  35 #include "AliESDEvent.h"
  36 #include "AliESDtrack.h"
  37 #include "AliESDVertex.h"
  38 #include "AliESDv0.h"
  39 #include "AliHLTMessage.h"
  40 #include "TMath.h"
  41 #include "AliKFParticle.h"
  42 #include "AliKFVertex.h"
  43
  44
  45 /** ROOT macro for the implementation of ROOT specific class methods */
  46 ClassImp(AliHLTGlobalVertexerComponent)
  47
  48 AliHLTGlobalVertexerComponent::AliHLTGlobalVertexerComponent()
  49 :
  50   fESD(0),
  51   fTrackInfos(0),
  52   fPrimaryVtx(),
  53   fHistPrimVertexXY(0),
  54   fHistPrimVertexZX(0),
  55   fHistPrimVertexZY(0),
  56   fNEvents(0),
  57   fPlotHistograms(1),
  58   fFitTracksToVertex(1),
  59   fConstrainedTrackDeviation(4.),
  60   fV0DaughterPrimDeviation( 2.5 ),
  61   fV0PrimDeviation( 3.5 ),
  62   fV0Chi(3.5),
  63   fV0DecayLengthInSigmas(3.)
  64 {
  65   // see header file for class documentation
  66   // or
  67   // refer to README to build package
  68   // or
  69   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  70
  71 }
  72
  73 AliHLTGlobalVertexerComponent::~AliHLTGlobalVertexerComponent()
  74 {
  75   // see header file for class documentation
  76
  77   if( fTrackInfos ) delete[] fTrackInfos;
  78   if( fHistPrimVertexXY ) delete   fHistPrimVertexXY;
  79   if( fHistPrimVertexZX ) delete   fHistPrimVertexZX;
  80   if( fHistPrimVertexZY ) delete   fHistPrimVertexZY;
  81 }
  82
  83 // Public functions to implement AliHLTComponent's interface.
  84 // These functions are required for the registration process
  85
  86 const char* AliHLTGlobalVertexerComponent::GetComponentID()
  87 {
  88   // see header file for class documentation
  89
  90   return "GlobalVertexer";
  91 }
  92
  93 void AliHLTGlobalVertexerComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
  94 {
  95   // see header file for class documentation
  96   list.clear();
  97   list.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut );
  98 }
  99
 100 AliHLTComponentDataType AliHLTGlobalVertexerComponent::GetOutputDataType()
 101 {
 102   // see header file for class documentation
 103   return kAliHLTMultipleDataType;
 104 }
 105
 106 int AliHLTGlobalVertexerComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
 107
 108 {
 109   // see header file for class documentation
 110   tgtList.clear();
 111   tgtList.push_back(kAliHLTDataTypeESDObject|kAliHLTDataOriginOut);
 112   tgtList.push_back(kAliHLTDataTypeHistogram);
 113   return tgtList.size();
 114 }
 115
 116 void AliHLTGlobalVertexerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
 117 {
 118   // see header file for class documentation
 119   // XXX TODO: Find more realistic values.
 120   constBase = 80000;
 121   inputMultiplier = 2.;
 122 }
 123
 124 AliHLTComponent* AliHLTGlobalVertexerComponent::Spawn()
 125 {
 126   // see header file for class documentation
 127   return new AliHLTGlobalVertexerComponent;
 128 }
 129
 130 int AliHLTGlobalVertexerComponent::DoInit( int argc, const char** argv )
 131 {
 132   // init
 133
 134   fHistPrimVertexXY = new TH2F("primVertexXY","HLT:  Primary vertex distribution in XY",60,-2.,2.,60,-2.,2.);
 135   fHistPrimVertexXY->SetMarkerStyle(8);
 136   fHistPrimVertexXY->SetMarkerSize(0.4);
 137   fHistPrimVertexXY->SetYTitle("Y");
 138   fHistPrimVertexXY->SetXTitle("X");
 139   fHistPrimVertexXY->SetStats(0);
 140   //fHistPrimVertexXY->SetOption("COLZ");
 141
 142   fHistPrimVertexZX = new TH2F("primVertexZX","HLT:  Primary vertex distribution in ZX",60,-15.,15.,60,-2.,2.);
 143   fHistPrimVertexZX->SetMarkerStyle(8);
 144   fHistPrimVertexZX->SetMarkerSize(0.4);
 145   fHistPrimVertexZX->SetYTitle("X");
 146   fHistPrimVertexZX->SetXTitle("Z");
 147   fHistPrimVertexZX->SetStats(0);
 148   //fHistPrimVertexZX->SetOption("COLZ");
 149
 150   fHistPrimVertexZY = new TH2F("primVertexZY","HLT:  Primary vertex distribution in ZY",60,-10.,10.,60,-2.,2.);
 151   fHistPrimVertexZY->SetMarkerStyle(8);
 152   fHistPrimVertexZY->SetMarkerSize(0.4);
 153   fHistPrimVertexZY->SetYTitle("Y");
 154   fHistPrimVertexZY->SetXTitle("Z");
 155   fHistPrimVertexZY->SetStats(0);
 156   //fHistPrimVertexZY->SetOption("COLZ");
 157
 158   fNEvents =0;
 159
 160   int iResult=0;
 161   TString configuration="";
 162   TString argument="";
 163   for (int i=0; i<argc && iResult>=0; i++) {
 164     argument=argv[i];
 165     if (!configuration.IsNull()) configuration+=" ";
 166     configuration+=argument;
 167   }
 168
 169   if (!configuration.IsNull()) {
 170     iResult=Configure(configuration.Data());
 171   }
 172
 173   return iResult;
 174 }
 175
 176 int AliHLTGlobalVertexerComponent::DoDeinit()
 177 {
 178   // see header file for class documentation
 179
 180   if( fTrackInfos ) delete[] fTrackInfos;
 181   if( fHistPrimVertexXY ) delete fHistPrimVertexXY;
 182   if( fHistPrimVertexZX ) delete fHistPrimVertexZX;
 183   if( fHistPrimVertexZY ) delete fHistPrimVertexZY;
 184
 185   fTrackInfos = 0;
 186   fHistPrimVertexXY = 0;
 187   fHistPrimVertexZX = 0;
 188   fHistPrimVertexZY = 0;
 189
 190   fPlotHistograms = 1;
 191   fFitTracksToVertex = 1;
 192   fConstrainedTrackDeviation = 4.;
 193   fV0DaughterPrimDeviation = 2.5 ;
 194   fV0PrimDeviation =3.5;
 195   fV0Chi = 3.5;
 196   fV0DecayLengthInSigmas = 3.;
 197   fNEvents = 0;
 198   return 0;
 199 }
 200
 201 int AliHLTGlobalVertexerComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/)
 202 {
 203
 204   if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
 205     return 0;
 206
 207   int iResult = 0;
 208
 209   fNEvents++;
 210
 211   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
 212
 213     AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
 214     if( !event ) continue;
 215     event->GetStdContent();
 216
 217     // primary vertex & V0's
 218
 219     SetESD( event );
 220     FindPrimaryVertex();
 221     FindV0s();
 222     const AliESDVertex *vPrim = event->GetPrimaryVertexTracks();
 223
 224     // plot histograms
 225
 226     if( fPlotHistograms ){
 227       if( vPrim && vPrim->GetNContributors()>=3 ){
 228         if( fHistPrimVertexXY ) fHistPrimVertexXY->Fill( vPrim->GetX(),vPrim->GetY()  );
 229         if( fHistPrimVertexZX ) fHistPrimVertexZX->Fill( vPrim->GetZ(),vPrim->GetX()  );
 230         if( fHistPrimVertexZY ) fHistPrimVertexZY->Fill( vPrim->GetZ(),vPrim->GetY()  );
 231       }
 232       if( fHistPrimVertexXY ) PushBack( fHistPrimVertexXY, kAliHLTDataTypeHistogram,0);
 233       if( fHistPrimVertexZX ) PushBack( fHistPrimVertexZX, kAliHLTDataTypeHistogram,0);
 234       if( fHistPrimVertexZY ) PushBack( fHistPrimVertexZY, kAliHLTDataTypeHistogram,0);
 235     }
 236
 237     iResult = PushBack( event, kAliHLTDataTypeESDObject|kAliHLTDataOriginOut, 0);
 238     if( iResult<0 ) break;
 239   }
 240
 241   return iResult;
 242 }
 243
 244 int AliHLTGlobalVertexerComponent::Configure(const char* arguments)
 245 {
 246
 247   int iResult=0;
 248   if (!arguments) return iResult;
 249
 250   TString allArgs=arguments;
 251   TString argument;
 252
 253   TObjArray* pTokens=allArgs.Tokenize(" ");
 254   int bMissingParam=0;
 255
 256   if (pTokens) {
 257     for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
 258       argument=((TObjString*)pTokens->At(i))->GetString();
 259       if (argument.IsNull()) continue;
 260
 261       if (argument.CompareTo("-fitTracksToVertex")==0) {
 262         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 263         HLTInfo("fitTracksToVertex is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 264         fFitTracksToVertex=((TObjString*)pTokens->At(i))->GetString().Atoi();
 265         continue;
 266       }
 267       else if (argument.CompareTo("-plotHistograms")==0) {
 268         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 269         HLTInfo("plotHistograms is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 270         fPlotHistograms=((TObjString*)pTokens->At(i))->GetString().Atoi();
 271         continue;
 272       }
 273       else if (argument.CompareTo("-constrainedTrackDeviation")==0) {
 274         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 275         HLTInfo("constrainedTrackDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 276         fConstrainedTrackDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 277         continue;
 278       }
 279       else if (argument.CompareTo("-v0DaughterPrimDeviation")==0) {
 280         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 281         HLTInfo("v0DaughterPrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 282         fV0DaughterPrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 283         continue;
 284       }
 285       else if (argument.CompareTo("-v0PrimDeviation")==0) {
 286         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 287         HLTInfo("v0PrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 288         fV0PrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
 289         continue;
 290       }
 291       else if (argument.CompareTo("-v0Chi")==0) {
 292         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 293         HLTInfo("v0Chi is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 294         fV0Chi=((TObjString*)pTokens->At(i))->GetString().Atof();
 295         continue;
 296       }
 297       else if (argument.CompareTo("-v0DecayLengthInSigmas")==0) {
 298         if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
 299         HLTInfo("v0DecayLengthInSigmas is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
 300         fV0DecayLengthInSigmas=((TObjString*)pTokens->At(i))->GetString().Atof();
 301         continue;
 302       }
 303       else {
 304         HLTError("unknown argument %s", argument.Data());
 305         iResult=-EINVAL;
 306         break;
 307       }
 308     }
 309     delete pTokens;
 310   }
 311   if (bMissingParam) {
 312     HLTError("missing parameter for argument %s", argument.Data());
 313     iResult=-EINVAL;
 314   }
 315
 316   return iResult;
 317 }
 318
 319 int AliHLTGlobalVertexerComponent::Reconfigure(const char* cdbEntry, const char* chainId)
 320 {
 321   // see header file for class documentation
 322   int iResult=0;
 323   const char* path="HLT/ConfigTPC/KryptonHistoComponent";
 324   const char* defaultNotify="";
 325   if (cdbEntry) {
 326     path=cdbEntry;
 327     defaultNotify=" (default)";
 328   }
 329   if (path) {
 330     HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
 331     AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path/*,GetRunNo()*/);
 332     if (pEntry) {
 333       TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
 334       if (pString) {
 335         HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
 336         iResult=Configure(pString->GetString().Data());
 337       } else {
 338         HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
 339       }
 340     } else {
 341       HLTError("can not fetch object \"%s\" from CDB", path);
 342     }
 343   }
 344
 345   return iResult;
 346 }
 347
 348
 349
 350 void AliHLTGlobalVertexerComponent::SetESD( AliESDEvent *event )
 351 {
 352   //* Fill fTrackInfo array
 353
 354   if( fTrackInfos ) delete[] fTrackInfos;
 355   fTrackInfos = 0;
 356   fESD = event;
 357
 358   AliKFParticle::SetField( fESD->GetMagneticField() );
 359
 360   Int_t nESDTracks=event->GetNumberOfTracks();
 361
 362   fTrackInfos = new AliESDTrackInfo[ nESDTracks ];
 363
 364   for (Int_t iTr=0; iTr<nESDTracks; iTr++){
 365
 366     AliESDTrackInfo &info = fTrackInfos[iTr];
 367     info.fOK = 0;
 368     info.fPrimUsedFlag = 0;
 369
 370     //* track quality check
 371
 372     AliESDtrack *pTrack = event->GetTrack(iTr);
 373     if( !pTrack  ) continue;
 374     if (pTrack->GetKinkIndex(0)>0) continue;
 375     if ( !( pTrack->GetStatus()&AliESDtrack::kTPCin ) ) continue;
 376
 377     //* Construct KFParticle for the track
 378
 379     if(  pTrack->GetStatus()&AliESDtrack::kITSin ){
 380       info.fParticle = AliKFParticle( *pTrack, 211 );
 381     } else {
 382       info.fParticle = AliKFParticle( *pTrack->GetInnerParam(), 211 );
 383     }
 384     info.fOK = 1;
 385   }
 386 }
 387
 388
 389 void AliHLTGlobalVertexerComponent::FindPrimaryVertex(  )
 390 {
 391   //* Find event primary vertex
 392
 393   int nTracks = fESD->GetNumberOfTracks();
 394
 395   const AliKFParticle **vSelected = new const AliKFParticle*[nTracks]; //* Selected particles for vertex fit
 396   Int_t *vIndex = new int [nTracks];                    //* Indices of selected particles
 397   Bool_t *vFlag = new bool [nTracks];                    //* Flags returned by the vertex finder
 398
 399   fPrimaryVtx.Initialize();
 400   fPrimaryVtx.SetBeamConstraint(fESD->GetDiamondX(),fESD->GetDiamondY(),0,
 401                                 TMath::Sqrt(fESD->GetSigma2DiamondX()),TMath::Sqrt(fESD->GetSigma2DiamondY()),5.3);
 402
 403   Int_t nSelected = 0;
 404   for( Int_t i = 0; i<nTracks; i++){
 405     if(!fTrackInfos[i].fOK ) continue;
 406     //if( fESD->GetTrack(i)->GetTPCNcls()<60  ) continue;
 407     const AliKFParticle &p = fTrackInfos[i].fParticle;
 408     Double_t chi = p.GetDeviationFromVertex( fPrimaryVtx );
 409     if( chi > fConstrainedTrackDeviation ) continue;
 410     vSelected[nSelected] = &(fTrackInfos[i].fParticle);
 411     vIndex[nSelected] = i;
 412     nSelected++;
 413   }
 414   fPrimaryVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, fConstrainedTrackDeviation );
 415   for( Int_t i = 0; i<nSelected; i++){
 416     if( vFlag[i] ) fTrackInfos[vIndex[i]].fPrimUsedFlag = 1;
 417   }
 418
 419   for( Int_t i = 0; i<nTracks; i++ ){
 420     AliESDTrackInfo &info = fTrackInfos[i];
 421     info.fPrimDeviation = info.fParticle.GetDeviationFromVertex( fPrimaryVtx );
 422   }
 423
 424   if( fPrimaryVtx.GetNContributors()>=3 ){
 425     AliESDVertex vESD( fPrimaryVtx.Parameters(), fPrimaryVtx.CovarianceMatrix(), fPrimaryVtx.GetChi2(), fPrimaryVtx.GetNContributors() );
 426     fESD->SetPrimaryVertexTracks( &vESD );
 427
 428     // relate the tracks to vertex
 429
 430     if( fFitTracksToVertex ){
 431       for( Int_t i = 0; i<nTracks; i++ ){
 432         if( !fTrackInfos[i].fPrimUsedFlag ) continue;
 433         if( fTrackInfos[i].fPrimDeviation > fConstrainedTrackDeviation ) continue;
 434         fESD->GetTrack(i)->RelateToVertex( &vESD, fESD->GetMagneticField(),100. );
 435       }
 436     }
 437
 438   } else {
 439     for( Int_t i = 0; i<nTracks; i++)
 440       fTrackInfos[i].fPrimUsedFlag = 0;
 441   }
 442
 443
 444   delete[] vSelected;
 445   delete[] vIndex;
 446   delete[] vFlag;
 447 }
 448
 449
 450 void AliHLTGlobalVertexerComponent::FindV0s(  )
 451 {
 452   //* V0 finder
 453
 454   int nTracks = fESD->GetNumberOfTracks();
 455   //AliKFVertex primVtx( *fESD->GetPrimaryVertexTracks() );
 456   AliKFVertex &primVtx = fPrimaryVtx;
 457   if( primVtx.GetNContributors()<3 ) return;
 458
 459   bool *constrainedV0   = new bool[nTracks];
 460   for( Int_t iTr = 0; iTr<nTracks; iTr++ ){
 461     constrainedV0[iTr] = 0;
 462   }
 463
 464
 465   for( Int_t iTr = 0; iTr<nTracks; iTr++ ){ //* first daughter
 466
 467     AliESDTrackInfo &info = fTrackInfos[iTr];
 468     if( !info.fOK ) continue;
 469     if( info.fParticle.GetQ() >0 ) continue;
 470     if( info.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
 471
 472     for( Int_t jTr = 0; jTr<nTracks; jTr++ ){  //* second daughter
 473
 474       AliESDTrackInfo &jnfo = fTrackInfos[jTr];
 475       if( !jnfo.fOK ) continue;
 476       if( jnfo.fParticle.GetQ() < 0 ) continue;
 477       if( jnfo.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
 478
 479       //* construct V0 mother
 480
 481       AliKFParticle v0( info.fParticle, jnfo.fParticle );
 482
 483       //* check V0 Chi^2
 484
 485       if( v0.GetChi2()<0 || v0.GetChi2() > fV0Chi*fV0Chi*v0.GetNDF() ) continue;
 486
 487       //* subtruct daughters from primary vertex
 488
 489       AliKFVertex primVtxCopy = primVtx;
 490
 491       if( info.fPrimUsedFlag ){
 492         if( primVtxCopy.GetNContributors()<=2 ) continue;
 493         primVtxCopy -= info.fParticle;
 494       }
 495       if( jnfo.fPrimUsedFlag ){
 496         if( primVtxCopy.GetNContributors()<=2 ) continue;
 497         primVtxCopy -= jnfo.fParticle;
 498       }
 499       //* Check v0 Chi^2 deviation from primary vertex
 500
 501       if( v0.GetDeviationFromVertex( primVtxCopy ) > fV0PrimDeviation ) continue;
 502
 503       //* Add V0 to primary vertex to improve the primary vertex resolution
 504
 505       primVtxCopy += v0;
 506
 507       //* Set production vertex for V0
 508
 509       v0.SetProductionVertex( primVtxCopy );
 510
 511       //* Get V0 decay length with estimated error
 512
 513       Double_t length, sigmaLength;
 514       if( v0.GetDecayLength( length, sigmaLength ) ) continue;
 515
 516       //* Reject V0 if it decays too close[sigma] to the primary vertex
 517
 518       if( length  < fV0DecayLengthInSigmas*sigmaLength ) continue;
 519
 520       //* add ESD v0
 521
 522       AliESDv0 v0ESD( *fESD->GetTrack( iTr ), iTr, *fESD->GetTrack( jTr ), jTr );
 523       fESD->AddV0( &v0ESD );
 524
 525       // relate the tracks to vertex
 526
 527       if( fFitTracksToVertex ){
 528         if( constrainedV0[iTr] || constrainedV0[jTr]
 529             || info.fPrimDeviation < fConstrainedTrackDeviation || jnfo.fPrimDeviation < fConstrainedTrackDeviation ) continue;
 530         AliESDVertex vESD(v0.Parameters(), v0.CovarianceMatrix(), v0.GetChi2(), 2);
 531         fESD->GetTrack(iTr)->RelateToVertex( &vESD, fESD->GetMagneticField(),100. );
 532         fESD->GetTrack(jTr)->RelateToVertex( &vESD, fESD->GetMagneticField(),100. );
 533         constrainedV0[iTr] = 1;
 534         constrainedV0[jTr] = 1;
 535       }
 536     }
 537   }
 538   delete[] constrainedV0;
 539 }
 540