[u/mrichter/AliRoot.git] / HLT / global / AliHLTGlobalVertexerComponent.cxx

//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: S.Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>    *
//*                  for The ALICE HLT Project.                            *
//*                                                                        *
//* Permission to use, copy, modify and distribute this software and its   *
//* documentation strictly for non-commercial purposes is hereby granted   *
//* without fee, provided that the above copyright notice appears in all   *
//* copies and that both the copyright notice and this permission notice   *
//* appear in the supporting documentation. The authors make no claims     *
//* about the suitability of this software for any purpose. It is          *
//* provided "as is" without express or implied warranty.                  *
//**************************************************************************

/** @file   AliHLTGlobalVertexerComponent.cxx
    @author Sergey Gorbunov
    @brief  Component for reconstruct primary vertex and V0's
*/

#if __GNUC__>= 3
using namespace std;
#endif

#include "AliHLTGlobalVertexerComponent.h"
#include "AliCDBEntry.h"
#include "AliCDBManager.h"
#include <TFile.h>
#include <TString.h>
#include "TObjString.h"
#include "TObjArray.h"
#include "TH1F.h"
#include "TH2F.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliESDv0.h"
#include "AliHLTMessage.h"
#include "TMath.h"
#include "AliKFParticle.h"
#include "AliKFVertex.h"


/** ROOT macro for the implementation of ROOT specific class methods */
ClassImp(AliHLTGlobalVertexerComponent)

AliHLTGlobalVertexerComponent::AliHLTGlobalVertexerComponent()
:
  fESD(0),
  fTrackInfos(0),
  fPrimaryVtx(),  
  fHistPrimVertexXY(0),
  fHistPrimVertexZX(0),
  fHistPrimVertexZY(0),
  fNEvents(0),
  fPlotHistograms(0),
  fFitTracksToVertex(1),
  fConstrainedTrackDeviation(4.),
  fV0DaughterPrimDeviation( 2.5 ),
  fV0PrimDeviation( 3.5 ),
  fV0Chi(3.5),
  fV0DecayLengthInSigmas(3.)
{
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt

}

AliHLTGlobalVertexerComponent::~AliHLTGlobalVertexerComponent()
{
  // see header file for class documentation

  if( fTrackInfos ) delete[] fTrackInfos;
  if( fHistPrimVertexXY ) delete   fHistPrimVertexXY;
  if( fHistPrimVertexZX ) delete   fHistPrimVertexZX;
  if( fHistPrimVertexZY ) delete   fHistPrimVertexZY;
}

// Public functions to implement AliHLTComponent's interface.
// These functions are required for the registration process

const char* AliHLTGlobalVertexerComponent::GetComponentID()
{
  // see header file for class documentation
  
  return "GlobalVertexer";
}

void AliHLTGlobalVertexerComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
{
  // see header file for class documentation
  list.clear();
  list.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut );
}

AliHLTComponentDataType AliHLTGlobalVertexerComponent::GetOutputDataType()
{
  // see header file for class documentation
  return kAliHLTMultipleDataType;
}

int AliHLTGlobalVertexerComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)

{
  // see header file for class documentation
  tgtList.clear();
  tgtList.push_back(kAliHLTDataTypeESDObject|kAliHLTDataOriginOut);
  tgtList.push_back(kAliHLTDataTypeHistogram);
  return tgtList.size();
}

void AliHLTGlobalVertexerComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
{
  // see header file for class documentation
  // XXX TODO: Find more realistic values.
  constBase = 80000;
  inputMultiplier = 2.;
}

AliHLTComponent* AliHLTGlobalVertexerComponent::Spawn()
{
  // see header file for class documentation
  return new AliHLTGlobalVertexerComponent;
}

int AliHLTGlobalVertexerComponent::DoInit( int argc, const char** argv )
{
  // init
  
  fHistPrimVertexXY = new TH2F("primVertexXY","HLT:  Primary vertex distribution in XY",60,-2.,2.,60,-2.,2.);
  fHistPrimVertexXY->SetMarkerStyle(8);
  fHistPrimVertexXY->SetMarkerSize(0.4);
  fHistPrimVertexXY->SetYTitle("Y");
  fHistPrimVertexXY->SetXTitle("X");
  fHistPrimVertexXY->SetStats(0);
  //fHistPrimVertexXY->SetOption("COLZ");

  fHistPrimVertexZX = new TH2F("primVertexZX","HLT:  Primary vertex distribution in ZX",60,-15.,15.,60,-2.,2.);
  fHistPrimVertexZX->SetMarkerStyle(8);
  fHistPrimVertexZX->SetMarkerSize(0.4);
  fHistPrimVertexZX->SetYTitle("X");
  fHistPrimVertexZX->SetXTitle("Z");
  fHistPrimVertexZX->SetStats(0);
  //fHistPrimVertexZX->SetOption("COLZ");

  fHistPrimVertexZY = new TH2F("primVertexZY","HLT:  Primary vertex distribution in ZY",60,-10.,10.,60,-2.,2.);
  fHistPrimVertexZY->SetMarkerStyle(8);
  fHistPrimVertexZY->SetMarkerSize(0.4);
  fHistPrimVertexZY->SetYTitle("Y");
  fHistPrimVertexZY->SetXTitle("Z");
  fHistPrimVertexZY->SetStats(0);
  //fHistPrimVertexZY->SetOption("COLZ");

  fNEvents =0;

  int iResult=0;
  TString configuration="";
  TString argument="";
  for (int i=0; i<argc && iResult>=0; i++) {
    argument=argv[i];
    if (!configuration.IsNull()) configuration+=" ";
    configuration+=argument;
  }
  
  if (!configuration.IsNull()) {
    iResult=Configure(configuration.Data());
  }  

  return iResult; 
}
  
int AliHLTGlobalVertexerComponent::DoDeinit()
{
  // see header file for class documentation

  if( fTrackInfos ) delete[] fTrackInfos;
  if( fHistPrimVertexXY ) delete fHistPrimVertexXY;
  if( fHistPrimVertexZX ) delete fHistPrimVertexZX;
  if( fHistPrimVertexZY ) delete fHistPrimVertexZY;

  fTrackInfos = 0;
  fHistPrimVertexXY = 0;
  fHistPrimVertexZX = 0;
  fHistPrimVertexZY = 0;

  fPlotHistograms = 0;
  fFitTracksToVertex = 1;
  fConstrainedTrackDeviation = 4.;
  fV0DaughterPrimDeviation = 2.5 ;
  fV0PrimDeviation =3.5;
  fV0Chi = 3.5;
  fV0DecayLengthInSigmas = 3.;
  fNEvents = 0;
  return 0;
}

int AliHLTGlobalVertexerComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/, AliHLTComponentTriggerData& /*trigData*/)
{
  
  if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
    return 0;

  int iResult = 0;

  fNEvents++;

  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {

    AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
    if( !event ) continue;
    event->GetStdContent();

    // primary vertex & V0's 
          
    SetESD( event );
    FindPrimaryVertex();
    FindV0s();
    const AliESDVertex *vPrim = event->GetPrimaryVertexTracks();

    // plot histograms

    if( fPlotHistograms ){
      if( vPrim && vPrim->GetNContributors()>=3 ){
	if( fHistPrimVertexXY ) fHistPrimVertexXY->Fill( vPrim->GetX(),vPrim->GetY()  );
	if( fHistPrimVertexZX ) fHistPrimVertexZX->Fill( vPrim->GetZ(),vPrim->GetX()  );
	if( fHistPrimVertexZY ) fHistPrimVertexZY->Fill( vPrim->GetZ(),vPrim->GetY()  );    
      }  
      if( fHistPrimVertexXY ) PushBack( fHistPrimVertexXY, kAliHLTDataTypeHistogram,0);
      if( fHistPrimVertexZX ) PushBack( fHistPrimVertexZX, kAliHLTDataTypeHistogram,0);
      if( fHistPrimVertexZY ) PushBack( fHistPrimVertexZY, kAliHLTDataTypeHistogram,0);
    }
    
    iResult = PushBack( event, kAliHLTDataTypeESDObject|kAliHLTDataOriginOut, 0);
    if( iResult<0 ) break;
  }
    
  return iResult;
}

int AliHLTGlobalVertexerComponent::Configure(const char* arguments)
{
  
  int iResult=0;
  if (!arguments) return iResult;
  
  TString allArgs=arguments;
  TString argument;
  
  TObjArray* pTokens=allArgs.Tokenize(" ");
  int bMissingParam=0;

  if (pTokens) {
    for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
      argument=((TObjString*)pTokens->At(i))->GetString();
      if (argument.IsNull()) continue;
      
      if (argument.CompareTo("-fitTracksToVertex")==0) {
	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	HLTInfo("fitTracksToVertex is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
	fFitTracksToVertex=((TObjString*)pTokens->At(i))->GetString().Atoi();
	continue;
      }
      else if (argument.CompareTo("-plotHistograms")==0) {
	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	HLTInfo("plotHistograms is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
	fPlotHistograms=((TObjString*)pTokens->At(i))->GetString().Atoi();
	continue;
      }
      else if (argument.CompareTo("-constrainedTrackDeviation")==0) {
	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	HLTInfo("constrainedTrackDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
	fConstrainedTrackDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
	continue;
      }
      else if (argument.CompareTo("-v0DaughterPrimDeviation")==0) {
	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	HLTInfo("v0DaughterPrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
	fV0DaughterPrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
	continue;
      }
      else if (argument.CompareTo("-v0PrimDeviation")==0) {
	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	HLTInfo("v0PrimDeviation is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
	fV0PrimDeviation=((TObjString*)pTokens->At(i))->GetString().Atof();
	continue;
      }
      else if (argument.CompareTo("-v0Chi")==0) {
	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	HLTInfo("v0Chi is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
	fV0Chi=((TObjString*)pTokens->At(i))->GetString().Atof();
	continue;
      }
      else if (argument.CompareTo("-v0DecayLengthInSigmas")==0) {
	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	HLTInfo("v0DecayLengthInSigmas is set set to: %s", ((TObjString*)pTokens->At(i))->GetString().Data());
	fV0DecayLengthInSigmas=((TObjString*)pTokens->At(i))->GetString().Atof();
	continue;
      }
      else {
	HLTError("unknown argument %s", argument.Data());
	iResult=-EINVAL;
	break;
      }
    }
    delete pTokens;
  }
  if (bMissingParam) {
    HLTError("missing parameter for argument %s", argument.Data());
    iResult=-EINVAL;
  }  
  
  return iResult;
}

int AliHLTGlobalVertexerComponent::Reconfigure(const char* cdbEntry, const char* chainId)
{
  // see header file for class documentation
  int iResult=0;
  const char* path="HLT/ConfigTPC/KryptonHistoComponent";
  const char* defaultNotify="";
  if (cdbEntry) {
    path=cdbEntry;
    defaultNotify=" (default)";
  }
  if (path) {
    HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
    AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path/*,GetRunNo()*/);
    if (pEntry) {
      TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
      if (pString) {
	HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
	iResult=Configure(pString->GetString().Data());
      } else {
	HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
      }
    } else {
      HLTError("can not fetch object \"%s\" from CDB", path);
    }
  }

  return iResult;
}


void AliHLTGlobalVertexerComponent::SetESD( AliESDEvent *event )
{
  //* Fill fTrackInfo array

  if( fTrackInfos ) delete[] fTrackInfos;
  fTrackInfos = 0;
  fESD = event;

  AliKFParticle::SetField( fESD->GetMagneticField() );

  Int_t nESDTracks=event->GetNumberOfTracks(); 

  fTrackInfos = new AliESDTrackInfo[ nESDTracks ];

  for (Int_t iTr=0; iTr<nESDTracks; iTr++){ 
  
    AliESDTrackInfo &info = fTrackInfos[iTr];
    info.fOK = 0;
    info.fPrimUsedFlag = 0;
    
    //* track quality check

    AliESDtrack *pTrack = event->GetTrack(iTr);    
    if( !pTrack  ) continue;
    if (pTrack->GetKinkIndex(0)>0) continue;
    if ( !( pTrack->GetStatus()&AliESDtrack::kTPCin ) ) continue;
    
    //* Construct KFParticle for the track

    if(  pTrack->GetStatus()&AliESDtrack::kITSin ){
      info.fParticle = AliKFParticle( *pTrack, 211 );    
    } else {
      info.fParticle = AliKFParticle( *pTrack->GetInnerParam(), 211 );    
    }
    info.fOK = 1;
  }
}


void AliHLTGlobalVertexerComponent::FindPrimaryVertex(  )
{
  //* Find event primary vertex

  int nTracks = fESD->GetNumberOfTracks();

  const AliKFParticle **vSelected = new const AliKFParticle*[nTracks]; //* Selected particles for vertex fit
  Int_t *vIndex = new int [nTracks];                    //* Indices of selected particles
  Bool_t *vFlag = new bool [nTracks];                    //* Flags returned by the vertex finder

  fPrimaryVtx.Initialize();
  fPrimaryVtx.SetBeamConstraint(fESD->GetDiamondX(),fESD->GetDiamondY(),0,
				TMath::Sqrt(fESD->GetSigma2DiamondX()),TMath::Sqrt(fESD->GetSigma2DiamondY()),5.3);
  
  Int_t nSelected = 0;
  for( Int_t i = 0; i<nTracks; i++){ 
    if(!fTrackInfos[i].fOK ) continue;
    //if( fESD->GetTrack(i)->GetTPCNcls()<60  ) continue;
    const AliKFParticle &p = fTrackInfos[i].fParticle;
    Double_t chi = p.GetDeviationFromVertex( fPrimaryVtx );      
    if( chi > fConstrainedTrackDeviation ) continue;
    vSelected[nSelected] = &(fTrackInfos[i].fParticle);
    vIndex[nSelected] = i;
    nSelected++;  
  }
  fPrimaryVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, fConstrainedTrackDeviation );
  for( Int_t i = 0; i<nSelected; i++){ 
    if( vFlag[i] ) fTrackInfos[vIndex[i]].fPrimUsedFlag = 1;
  }

  for( Int_t i = 0; i<nTracks; i++ ){
    AliESDTrackInfo &info = fTrackInfos[i];
    info.fPrimDeviation = info.fParticle.GetDeviationFromVertex( fPrimaryVtx );   
  }
  
  if( fPrimaryVtx.GetNContributors()>=3 ){
    AliESDVertex vESD( fPrimaryVtx.Parameters(), fPrimaryVtx.CovarianceMatrix(), fPrimaryVtx.GetChi2(), fPrimaryVtx.GetNContributors() );
    fESD->SetPrimaryVertexTracks( &vESD );

    // relate the tracks to vertex

    if( fFitTracksToVertex ){      
      for( Int_t i = 0; i<nTracks; i++ ){
	if( !fTrackInfos[i].fPrimUsedFlag ) continue;	  
	if( fTrackInfos[i].fPrimDeviation > fConstrainedTrackDeviation ) continue;
	fESD->GetTrack(i)->RelateToVertex( &vESD, fESD->GetMagneticField(),100. );
      }
    }

  } else {
    for( Int_t i = 0; i<nTracks; i++)
      fTrackInfos[i].fPrimUsedFlag = 0;
  }


  delete[] vSelected;
  delete[] vIndex;
  delete[] vFlag;
}


void AliHLTGlobalVertexerComponent::FindV0s(  )
{
  //* V0 finder

  int nTracks = fESD->GetNumberOfTracks();
  //AliKFVertex primVtx( *fESD->GetPrimaryVertexTracks() );
  AliKFVertex &primVtx = fPrimaryVtx;
  if( primVtx.GetNContributors()<3 ) return;

  bool *constrainedV0   = new bool[nTracks];
  for( Int_t iTr = 0; iTr<nTracks; iTr++ ){ 
    constrainedV0[iTr] = 0;
  }
  

  for( Int_t iTr = 0; iTr<nTracks; iTr++ ){ //* first daughter

    AliESDTrackInfo &info = fTrackInfos[iTr];
    if( !info.fOK ) continue;    
    if( info.fParticle.GetQ() >0 ) continue;    
    if( info.fPrimDeviation < fV0DaughterPrimDeviation ) continue;

    for( Int_t jTr = 0; jTr<nTracks; jTr++ ){  //* second daughter

      AliESDTrackInfo &jnfo = fTrackInfos[jTr];
      if( !jnfo.fOK ) continue;
      if( jnfo.fParticle.GetQ() < 0 ) continue;
      if( jnfo.fPrimDeviation < fV0DaughterPrimDeviation ) continue;
   
      //* construct V0 mother

      AliKFParticle v0( info.fParticle, jnfo.fParticle );     

      //* check V0 Chi^2
      
      if( v0.GetChi2()<0 || v0.GetChi2() > fV0Chi*fV0Chi*v0.GetNDF() ) continue;

      //* subtruct daughters from primary vertex 

      AliKFVertex primVtxCopy = primVtx;    
       
      if( info.fPrimUsedFlag ){	
	if( primVtxCopy.GetNContributors()<=2 ) continue;
	primVtxCopy -= info.fParticle;
      }
      if( jnfo.fPrimUsedFlag ){
	if( primVtxCopy.GetNContributors()<=2 ) continue;
	primVtxCopy -= jnfo.fParticle;
      }
      //* Check v0 Chi^2 deviation from primary vertex 

      if( v0.GetDeviationFromVertex( primVtxCopy ) > fV0PrimDeviation ) continue;

      //* Add V0 to primary vertex to improve the primary vertex resolution

      primVtxCopy += v0;      

      //* Set production vertex for V0

      v0.SetProductionVertex( primVtxCopy );

      //* Get V0 decay length with estimated error

      Double_t length, sigmaLength;
      if( v0.GetDecayLength( length, sigmaLength ) ) continue;

      //* Reject V0 if it decays too close[sigma] to the primary vertex

      if( length  < fV0DecayLengthInSigmas*sigmaLength ) continue;

      //* add ESD v0 
      
      AliESDv0 v0ESD( *fESD->GetTrack( iTr ), iTr, *fESD->GetTrack( jTr ), jTr );  
      fESD->AddV0( &v0ESD );

      // relate the tracks to vertex

      if( fFitTracksToVertex ){
	if( constrainedV0[iTr] || constrainedV0[jTr]
	    || info.fPrimDeviation < fConstrainedTrackDeviation || jnfo.fPrimDeviation < fConstrainedTrackDeviation ) continue;
	AliESDVertex vESD(v0.Parameters(), v0.CovarianceMatrix(), v0.GetChi2(), 2);
	fESD->GetTrack(iTr)->RelateToVertex( &vESD, fESD->GetMagneticField(),100. );
	fESD->GetTrack(jTr)->RelateToVertex( &vESD, fESD->GetMagneticField(),100. );
	constrainedV0[iTr] = 1;
	constrainedV0[jTr] = 1;	
      }
    }
  }
  delete[] constrainedV0;
}
Commit	Line	Data
4d5ee3db	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"
4d5ee3db	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),
5850d7e6	57	fPlotHistograms(0),
608cfbda	58	fFitTracksToVertex(1),
4d5ee3db	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
dc546924	77	if( fTrackInfos ) delete[] fTrackInfos;
	78	if( fHistPrimVertexXY ) delete fHistPrimVertexXY;
	79	if( fHistPrimVertexZX ) delete fHistPrimVertexZX;
	80	if( fHistPrimVertexZY ) delete fHistPrimVertexZY;
4d5ee3db	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
608cfbda	134	fHistPrimVertexXY = new TH2F("primVertexXY","HLT: Primary vertex distribution in XY",60,-2.,2.,60,-2.,2.);
4d5ee3db	135	fHistPrimVertexXY->SetMarkerStyle(8);
	136	fHistPrimVertexXY->SetMarkerSize(0.4);
	137	fHistPrimVertexXY->SetYTitle("Y");
	138	fHistPrimVertexXY->SetXTitle("X");
	139	fHistPrimVertexXY->SetStats(0);
608cfbda	140	//fHistPrimVertexXY->SetOption("COLZ");
4d5ee3db	141
608cfbda	142	fHistPrimVertexZX = new TH2F("primVertexZX","HLT: Primary vertex distribution in ZX",60,-15.,15.,60,-2.,2.);
4d5ee3db	143	fHistPrimVertexZX->SetMarkerStyle(8);
	144	fHistPrimVertexZX->SetMarkerSize(0.4);
	145	fHistPrimVertexZX->SetYTitle("X");
	146	fHistPrimVertexZX->SetXTitle("Z");
	147	fHistPrimVertexZX->SetStats(0);
608cfbda	148	//fHistPrimVertexZX->SetOption("COLZ");
4d5ee3db	149
608cfbda	150	fHistPrimVertexZY = new TH2F("primVertexZY","HLT: Primary vertex distribution in ZY",60,-10.,10.,60,-2.,2.);
4d5ee3db	151	fHistPrimVertexZY->SetMarkerStyle(8);
	152	fHistPrimVertexZY->SetMarkerSize(0.4);
	153	fHistPrimVertexZY->SetYTitle("Y");
	154	fHistPrimVertexZY->SetXTitle("Z");
	155	fHistPrimVertexZY->SetStats(0);
608cfbda	156	//fHistPrimVertexZY->SetOption("COLZ");
4d5ee3db	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
dc546924	180	if( fTrackInfos ) delete[] fTrackInfos;
	181	if( fHistPrimVertexXY ) delete fHistPrimVertexXY;
	182	if( fHistPrimVertexZX ) delete fHistPrimVertexZX;
	183	if( fHistPrimVertexZY ) delete fHistPrimVertexZY;
4d5ee3db	184
	185	fTrackInfos = 0;
	186	fHistPrimVertexXY = 0;
	187	fHistPrimVertexZX = 0;
	188	fHistPrimVertexZY = 0;
	189
5850d7e6	190	fPlotHistograms = 0;
608cfbda	191	fFitTracksToVertex = 1;
4d5ee3db	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
dc546924	207	int iResult = 0;
dc546924	208
4d5ee3db	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 ) );
dc546924	214	if( !event ) continue;
4d5ee3db	215	event->GetStdContent();
	216
	217	// primary vertex & V0's
608cfbda	218
	219	SetESD( event );
	220	FindPrimaryVertex();
	221	FindV0s();
4d5ee3db	222	const AliESDVertex *vPrim = event->GetPrimaryVertexTracks();
	223
	224	// plot histograms
	225
	226	if( fPlotHistograms ){
	227	if( vPrim && vPrim->GetNContributors()>=3 ){
dc546924	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() );
4d5ee3db	231	}
dc546924	232	if( fHistPrimVertexXY ) PushBack( fHistPrimVertexXY, kAliHLTDataTypeHistogram,0);
	233	if( fHistPrimVertexZX ) PushBack( fHistPrimVertexZX, kAliHLTDataTypeHistogram,0);
	234	if( fHistPrimVertexZY ) PushBack( fHistPrimVertexZY, kAliHLTDataTypeHistogram,0);
4d5ee3db	235	}
4d5ee3db	236
dc546924	237	iResult = PushBack( event, kAliHLTDataTypeESDObject\|kAliHLTDataOriginOut, 0);
dc546924	238	if( iResult<0 ) break;
4d5ee3db	239	}
4d5ee3db	240
dc546924	241	return iResult;
4d5ee3db	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());
608cfbda	264	fFitTracksToVertex=((TObjString*)pTokens->At(i))->GetString().Atoi();
4d5ee3db	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	}
4d5ee3db	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
dc546924	354	if( fTrackInfos ) delete[] fTrackInfos;
dc546924	355	fTrackInfos = 0;
4d5ee3db	356	fESD = event;
	357
	358	AliKFParticle::SetField( fESD->GetMagneticField() );
	359
	360	Int_t nESDTracks=event->GetNumberOfTracks();
dc546924	361
4d5ee3db	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
dc546924	379	if( pTrack->GetStatus()&AliESDtrack::kITSin ){
	380	info.fParticle = AliKFParticle( *pTrack, 211 );
	381	} else {
	382	info.fParticle = AliKFParticle( *pTrack->GetInnerParam(), 211 );
	383	}
4d5ee3db	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
608cfbda	430	if( fFitTracksToVertex ){
4d5ee3db	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() > fV0ChifV0Chiv0.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
608cfbda	527	if( fFitTracksToVertex ){
4d5ee3db	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