1 //**************************************************************************
2 //* This file is property of and copyright by the ALICE HLT Project *
3 //* ALICE Experiment at CERN, All rights reserved. *
5 //* Primary Authors: S.Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
6 //* for The ALICE HLT Project. *
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 //**************************************************************************
17 /** @file AliHLTGlobalVertexerComponent.cxx
18 @author Sergey Gorbunov
19 @brief Component for reconstruct primary vertex and V0's
26 #include "AliHLTGlobalVertexerComponent.h"
27 #include "AliHLTDataTypes.h"
28 #include "AliHLTExternalTrackParam.h"
29 #include "AliHLTGlobalBarrelTrack.h"
30 #include "AliCDBEntry.h"
31 #include "AliCDBManager.h"
34 #include "TObjString.h"
35 #include "TObjArray.h"
36 #include "AliESDEvent.h"
37 #include "AliESDtrack.h"
38 #include "AliESDVertex.h"
40 #include "AliHLTMessage.h"
42 #include "AliKFParticle.h"
43 #include "AliKFVertex.h"
44 #include "TStopwatch.h"
46 /** ROOT macro for the implementation of ROOT specific class methods */
47 ClassImp(AliHLTGlobalVertexerComponent)
49 AliHLTGlobalVertexerComponent::AliHLTGlobalVertexerComponent()
54 fFitTracksToVertex(1),
55 fConstrainedTrackDeviation(4.),
56 fV0DaughterPrimDeviation( 2.5 ),
57 fV0PrimDeviation( 3.5 ),
59 fV0DecayLengthInSigmas(3.),
61 fBenchmark("GlobalVertexer")
63 // see header file for class documentation
65 // refer to README to build package
67 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
71 AliHLTGlobalVertexerComponent::~AliHLTGlobalVertexerComponent()
73 // see header file for class documentation
75 if( fTrackInfos ) delete[] fTrackInfos;
78 // Public functions to implement AliHLTComponent's interface.
79 // These functions are required for the registration process
81 const char* AliHLTGlobalVertexerComponent::GetComponentID()
83 // see header file for class documentation
85 return "GlobalVertexer";
88 void AliHLTGlobalVertexerComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
90 // see header file for class documentation
92 list.push_back( kAliHLTDataTypeESDObject );
93 list.push_back( kAliHLTDataTypeTrack|kAliHLTDataOriginITS );
94 list.push_back( kAliHLTDataTypeTrack|kAliHLTDataOriginTPC );
97 AliHLTComponentDataType AliHLTGlobalVertexerComponent::GetOutputDataType()
99 // see header file for class documentation
100 return kAliHLTMultipleDataType;
103 int AliHLTGlobalVertexerComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
106 // see header file for class documentation
108 tgtList.push_back( kAliHLTDataTypeGlobalVertexer|kAliHLTDataOriginOut);
109 tgtList.push_back( kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut);
110 tgtList.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut);
111 return tgtList.size();
114 void AliHLTGlobalVertexerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
116 // see header file for class documentation
117 // XXX TODO: Find more realistic values.
119 inputMultiplier = 2.;
122 AliHLTComponent* AliHLTGlobalVertexerComponent::Spawn()
124 // see header file for class documentation
125 return new AliHLTGlobalVertexerComponent;
128 int AliHLTGlobalVertexerComponent::DoInit( int argc, const char** argv )
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;
140 TString configuration="";
142 for (int i=0; i<argc && iResult>=0; i++) {
144 if (!configuration.IsNull()) configuration+=" ";
145 configuration+=argument;
148 if (!configuration.IsNull()) {
149 iResult=Configure(configuration.Data());
155 int AliHLTGlobalVertexerComponent::DoDeinit()
157 // see header file for class documentation
159 if( fTrackInfos ) delete[] fTrackInfos;
162 fFitTracksToVertex = 1;
163 fConstrainedTrackDeviation = 4.;
164 fV0DaughterPrimDeviation = 2.5 ;
165 fV0PrimDeviation =3.5;
167 fV0DecayLengthInSigmas = 3.;
168 fV0TimeLimit = 10.e-3;
173 int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*evtData*/,
174 AliHLTComponentTriggerData& /*trigData*/ )
176 //cout<<"AliHLTGlobalVertexerComponent::DoEvent called"<<endl;
178 if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
182 //cout<<"SG: GlobalVertexer:: DoEvent called"<<endl;
183 fBenchmark.StartNewEvent();
186 vector< AliExternalTrackParam > tracks;
187 vector< int > trackId;
188 vector< pair<int,int> > v0s;
190 AliESDEvent *event = 0;
192 for( const AliHLTComponentBlockData *i= GetFirstInputBlock(kAliHLTDataTypeESDObject); i!=NULL; i=GetNextInputBlock() ){
193 fBenchmark.AddInput(i->fSize);
197 for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
198 event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
199 if( !event ) continue;
201 event->GetStdContent();
202 Int_t nESDTracks=event->GetNumberOfTracks();
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);
215 if( tracks.size()==0 ){
217 for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
218 pBlock!=NULL; pBlock=GetNextInputBlock()) {
220 fBenchmark.AddInput(pBlock->fSize);
222 AliHLTTracksData* dataPtr = reinterpret_cast<AliHLTTracksData*>( pBlock->fPtr );
223 int nTracks = dataPtr->fCount;
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 );
236 if( tracks.size()==0 ){
237 for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
238 pBlock!=NULL; pBlock=GetNextInputBlock()) {
240 fBenchmark.AddInput(pBlock->fSize);
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 );
256 // primary vertex & V0's
258 AliKFParticle::SetField( GetBz() );
262 { //* Fill fTrackInfo array
264 if( fTrackInfos ) delete[] fTrackInfos;
266 fNTracks=tracks.size();
267 fTrackInfos = new AliESDTrackInfo[ fNTracks ];
268 for (Int_t iTr=0; iTr<fNTracks; iTr++){
269 AliESDTrackInfo &info = fTrackInfos[iTr];
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;
286 int *buf = new int[sizeof(AliHLTGlobalVertexerData)/sizeof(int)+1 + fNTracks + 2*v0s.size()];
287 AliHLTGlobalVertexerData *data = reinterpret_cast<AliHLTGlobalVertexerData*>(buf);
289 if( data) { // fill the output structure
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];
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];
309 unsigned int blockSize = sizeof(AliHLTGlobalVertexerData) + (data->fNPrimTracks + 2*data->fNV0s)*sizeof(int);
311 iResult = PushBack( reinterpret_cast<void*>(data), blockSize, kAliHLTDataTypeGlobalVertexer|kAliHLTDataOriginOut );
312 fBenchmark.AddOutput(blockSize);
316 // output the vertex if found
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());
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());
335 HLTInfo(fBenchmark.GetStatistics());
340 int AliHLTGlobalVertexerComponent::Configure(const char* arguments)
344 if (!arguments) return iResult;
346 TString allArgs=arguments;
349 TObjArray* pTokens=allArgs.Tokenize(" ");
353 for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
354 argument=((TObjString*)pTokens->At(i))->GetString();
355 if (argument.IsNull()) continue;
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();
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();
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();
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();
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();
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();
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
401 HLTError("unknown argument %s", argument.Data());
409 HLTError("missing parameter for argument %s", argument.Data());
416 int AliHLTGlobalVertexerComponent::Reconfigure(const char* cdbEntry, const char* chainId)
418 // see header file for class documentation
420 return 0; // no CDB path is set so far
423 const char* path="HLT/ConfigTPC/KryptonHistoComponent";
424 const char* defaultNotify="";
427 defaultNotify=" (default)";
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/*,GetRunNo()*/);
433 TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
435 HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
436 iResult=Configure(pString->GetString().Data());
438 HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
441 HLTError("can not fetch object \"%s\" from CDB", path);
450 void AliHLTGlobalVertexerComponent::FindPrimaryVertex()
452 //* Find event primary vertex
454 const AliKFParticle **vSelected = new const AliKFParticle*[fNTracks]; //* Selected particles for the vertex fit
455 Int_t *vIndex = new int [fNTracks]; //* Indices of selected particles
456 Bool_t *vFlag = new bool [fNTracks]; //* Flags returned by the vertex finder
458 fPrimaryVtx.Initialize();
459 //fPrimaryVtx.SetBeamConstraint(fESD->GetDiamondX(),fESD->GetDiamondY(),0,
460 //TMath::Sqrt(fESD->GetSigma2DiamondX()),TMath::Sqrt(fESD->GetSigma2DiamondY()),5.3);
462 fPrimaryVtx.SetBeamConstraint( 0, 0, 0, 3., 3., 5.3 );
465 for( Int_t i = 0; i<fNTracks; i++){
466 if(!fTrackInfos[i].fOK ) continue;
467 //if( fESD->GetTrack(i)->GetTPCNcls()<60 ) continue;
468 const AliKFParticle &p = fTrackInfos[i].fParticle;
469 Double_t chi = p.GetDeviationFromVertex( fPrimaryVtx );
470 if( chi > fConstrainedTrackDeviation ) continue;
471 vSelected[nSelected] = &(fTrackInfos[i].fParticle);
472 vIndex[nSelected] = i;
476 fPrimaryVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, fConstrainedTrackDeviation );
478 for( Int_t i = 0; i<nSelected; i++){
479 if( vFlag[i] ) fTrackInfos[vIndex[i]].fPrimUsedFlag = 1;
482 for( Int_t i = 0; i<fNTracks; i++ ){
483 AliESDTrackInfo &info = fTrackInfos[i];
484 info.fPrimDeviation = info.fParticle.GetDeviationFromVertex( fPrimaryVtx );
488 if( fPrimaryVtx.GetNContributors()<3 ){
489 for( Int_t i = 0; i<fNTracks; i++)
490 fTrackInfos[i].fPrimUsedFlag = 0;
498 void AliHLTGlobalVertexerComponent::FindV0s( vector<pair<int,int> > &v0s )
502 AliKFVertex &primVtx = fPrimaryVtx;
503 if( primVtx.GetNContributors()<3 ) return;
509 for( Int_t iTr = 0; iTr<fNTracks && run; iTr++ ){ //* first daughter
511 AliESDTrackInfo &info = fTrackInfos[iTr];
512 if( !info.fOK ) continue;
513 if( info.fParticle.GetQ() >0 ) continue;
514 if( info.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
516 for( Int_t jTr = 0; jTr<fNTracks; jTr++ ){ //* second daughter
519 AliESDTrackInfo &jnfo = fTrackInfos[jTr];
520 if( !jnfo.fOK ) continue;
521 if( jnfo.fParticle.GetQ() < 0 ) continue;
522 if( jnfo.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
524 // check the time once a while...
526 if( (++statN)%100 ==0 ){
527 if( timer.RealTime()>= fV0TimeLimit ){ run = 0; break; }
531 //* check if the particles fit
533 if( info.fParticle.GetDeviationFromParticle(jnfo.fParticle) > fV0Chi ) continue;
535 //* construct V0 mother
537 AliKFParticle v0( info.fParticle, jnfo.fParticle );
541 if( v0.GetChi2()<0 || v0.GetChi2() > fV0Chi*fV0Chi*v0.GetNDF() ) continue;
543 //* subtruct daughters from primary vertex
545 AliKFVertex primVtxCopy = primVtx;
547 if( info.fPrimUsedFlag ){
548 if( primVtxCopy.GetNContributors()<=2 ) continue;
549 primVtxCopy -= info.fParticle;
551 if( jnfo.fPrimUsedFlag ){
552 if( primVtxCopy.GetNContributors()<=2 ) continue;
553 primVtxCopy -= jnfo.fParticle;
555 //* Check v0 Chi^2 deviation from primary vertex
557 if( v0.GetDeviationFromVertex( primVtxCopy ) > fV0PrimDeviation ) continue;
559 //* Add V0 to primary vertex to improve the primary vertex resolution
563 //* Set production vertex for V0
565 v0.SetProductionVertex( primVtxCopy );
567 //* Get V0 decay length with estimated error
569 Double_t length, sigmaLength;
570 if( v0.GetDecayLength( length, sigmaLength ) ) continue;
572 //* Reject V0 if it decays too close[sigma] to the primary vertex
574 if( length < fV0DecayLengthInSigmas*sigmaLength ) continue;
578 v0s.push_back(pair<int,int>(iTr,jTr));
586 void AliHLTGlobalVertexerComponent::FillESD( AliESDEvent *event, AliHLTGlobalVertexerData *data
589 //* put output of a vertexer to the esd event
591 Int_t nESDTracks = event->GetNumberOfTracks();
593 const int *listPrim = data->fTrackIndices;
594 const int *listV0 = data->fTrackIndices + data->fNPrimTracks;
596 std::map<int,int> mapId;
597 bool *constrainedToVtx = new bool[nESDTracks];
599 for( int i=0; i<nESDTracks; i++ ){
600 constrainedToVtx[i] = 0;
601 if( !event->GetTrack(i) ) continue;
602 mapId[ event->GetTrack(i)->GetID() ] = i;
605 if( data->fPrimNContributors >=3 ){
607 AliESDVertex vESD( data->fPrimP, data->fPrimC, data->fPrimChi2, data->fPrimNContributors );
608 event->SetPrimaryVertexTracks( &vESD );
610 // relate tracks to the primary vertex
612 if( data->fFitTracksToVertex ){
613 for( Int_t i = 0; i<data->fNPrimTracks; i++ ){
614 Int_t id = listPrim[ i ];
615 map<int,int>::iterator it = mapId.find(id);
616 if( it==mapId.end() ) continue;
617 Int_t itr = it->second;
618 event->GetTrack(itr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
619 constrainedToVtx[ itr ] = 1;
624 //* add ESD v0s and relate tracks to v0s
627 for( int i=0; i<data->fNV0s; i++ ){
629 Int_t id1 = listV0[ 2*i ];
630 Int_t id2 = listV0[ 2*i + 1];
631 map<int,int>::iterator it = mapId.find(id1);
632 if( it==mapId.end() ) continue;
633 Int_t iTr = it->second;
634 it = mapId.find(id2);
635 if( it==mapId.end() ) continue;
636 Int_t jTr = it->second;
638 AliESDv0 v0( *event->GetTrack( iTr ), iTr, *event->GetTrack( jTr ), jTr );
641 // relate the tracks to the vertex
643 if( data->fFitTracksToVertex ){
644 if( constrainedToVtx[iTr] || constrainedToVtx[jTr] ) continue;
646 double sigma[3] = {.1,.1,.1};
648 AliESDVertex vESD(pos, sigma);
649 event->GetTrack(iTr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
650 event->GetTrack(jTr)->RelateToVertex( &vESD, event->GetMagneticField(),100. );
651 constrainedToVtx[iTr] = 1;
652 constrainedToVtx[jTr] = 1;
656 delete[] constrainedToVtx;