]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/global/AliHLTGlobalVertexerComponent.cxx
Fix leaks in AliIsolationCut and AlianaOmegaToPi0Gamma
[u/mrichter/AliRoot.git] / HLT / global / AliHLTGlobalVertexerComponent.cxx
index 26027aa297dbbcc2ab63d9a3cb75406efa139cd1..0b86dba4f98d00f0d33182551d0ac90da424560d 100644 (file)
@@ -50,8 +50,7 @@ AliHLTGlobalVertexerComponent::AliHLTGlobalVertexerComponent()
 :
   fNTracks(0),
   fTrackInfos(0),
-  fPrimaryVtx(),  
-  fNEvents(0),
+  fPrimaryVtx(),   
   fFitTracksToVertex(1),
   fConstrainedTrackDeviation(4.),
   fV0DaughterPrimDeviation( 2.5 ),
@@ -59,15 +58,7 @@ AliHLTGlobalVertexerComponent::AliHLTGlobalVertexerComponent()
   fV0Chi(3.5),
   fV0DecayLengthInSigmas(3.),
   fV0TimeLimit(10.e-3), 
-  fStatTimeR( 0 ),
-  fStatTimeC( 0 ),
-  fStatTimeR1( 0 ),
-  fStatTimeC1( 0 ),
-  fStatTimeR2( 0 ),
-  fStatTimeC2( 0 ),
-  fStatTimeR3( 0 ),
-  fStatTimeC3( 0 ),
-  fStatNEvents(0)
+  fBenchmark("GlobalVertexer")
 {
   // see header file for class documentation
   // or
@@ -138,19 +129,13 @@ int AliHLTGlobalVertexerComponent::DoInit( int argc, const char** argv )
 {
   // init
 
-  fStatTimeR = 0;
-  fStatTimeC = 0;
-  fStatTimeR1 = 0;
-  fStatTimeC1 = 0;
-  fStatTimeR2 = 0;
-  fStatTimeC2 = 0;
-  fStatTimeR3 = 0;
-  fStatTimeC3 = 0;
-  fStatNEvents = 0;
+  fBenchmark.Reset();
+  fBenchmark.SetTimer(0,"total");
+  fBenchmark.SetTimer(1,"convert");
+  fBenchmark.SetTimer(2,"vprim");
+  fBenchmark.SetTimer(3,"v0");
   fV0TimeLimit = 10.e-3;
 
-  fNEvents =0;
-
   int iResult=0;
   TString configuration="";
   TString argument="";
@@ -180,7 +165,6 @@ int AliHLTGlobalVertexerComponent::DoDeinit()
   fV0PrimDeviation =3.5;
   fV0Chi = 3.5;
   fV0DecayLengthInSigmas = 3.;
-  fNEvents = 0;
   fV0TimeLimit = 10.e-3;
 
   return 0;
@@ -195,18 +179,25 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
     return 0;
 
   int iResult = 0;
-
-  fNEvents++;
-  TStopwatch timer;
+  //cout<<"SG: GlobalVertexer:: DoEvent called"<<endl;
+  fBenchmark.StartNewEvent();
+  fBenchmark.Start(0);
+  
   vector< AliExternalTrackParam > tracks;
   vector< int > trackId;
   vector< pair<int,int> > v0s;
 
   AliESDEvent *event = 0; 
 
+  for( const AliHLTComponentBlockData *i= GetFirstInputBlock(kAliHLTDataTypeESDObject); i!=NULL; i=GetNextInputBlock() ){
+    fBenchmark.AddInput(i->fSize);
+  }
+
+
   for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
     event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
     if( !event ) continue;
+    
     event->GetStdContent();
     Int_t nESDTracks=event->GetNumberOfTracks(); 
 
@@ -225,7 +216,9 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
     
     for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
         pBlock!=NULL; pBlock=GetNextInputBlock()) {
-      
+
+      fBenchmark.AddInput(pBlock->fSize);
+
       AliHLTTracksData* dataPtr = reinterpret_cast<AliHLTTracksData*>( pBlock->fPtr );
       int nTracks = dataPtr->fCount;
 
@@ -243,7 +236,9 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
   if( tracks.size()==0 ){
     for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
         pBlock!=NULL; pBlock=GetNextInputBlock()) {
-      
+
+      fBenchmark.AddInput(pBlock->fSize);
+
       AliHLTTracksData* dataPtr = reinterpret_cast<AliHLTTracksData*>( pBlock->fPtr );
       int nTracks = dataPtr->fCount;      
       AliHLTExternalTrackParam* currOutTrack = dataPtr->fTracklets;
@@ -262,10 +257,10 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
           
   AliKFParticle::SetField( GetBz() );
 
+  fBenchmark.Start(1);
+
   {  //* Fill fTrackInfo array
 
-    TStopwatch timer1;
-    
     if( fTrackInfos ) delete[] fTrackInfos;
     fTrackInfos = 0;
     fNTracks=tracks.size(); 
@@ -275,22 +270,18 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
       info.fOK = 1;
       info.fPrimUsedFlag = 0;
       info.fParticle = AliKFParticle( tracks[iTr], 211 );
+      for( int i=0; i<8; i++ ) if( !finite(info.fParticle.GetParameter(i)) ) info.fOK = 0;
+      for( int i=0; i<36; i++ ) if( !finite(info.fParticle.GetCovariance(i)) ) info.fOK = 0;
     }
-    timer1.Stop();
-    fStatTimeR1+=timer1.RealTime();
-    fStatTimeC1+=timer1.CpuTime();
   }
-  
-  TStopwatch timer2;
+
+  fBenchmark.Stop(1);
+  fBenchmark.Start(2);
   FindPrimaryVertex();
-  timer2.Stop();
-  fStatTimeR2+=timer2.RealTime();
-  fStatTimeC2+=timer2.CpuTime();
-  TStopwatch timer3;
+  fBenchmark.Stop(2);
+  fBenchmark.Start(3);
   FindV0s( v0s );
-  timer3.Stop();
-  fStatTimeR3+=timer3.RealTime();
-  fStatTimeC3+=timer3.CpuTime();
+  fBenchmark.Stop(3);
 
   int *buf = new int[sizeof(AliHLTGlobalVertexerData)/sizeof(int)+1 + fNTracks + 2*v0s.size()];
   AliHLTGlobalVertexerData *data = reinterpret_cast<AliHLTGlobalVertexerData*>(buf);
@@ -318,6 +309,7 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
     unsigned int blockSize = sizeof(AliHLTGlobalVertexerData) + (data->fNPrimTracks + 2*data->fNV0s)*sizeof(int);
 
     iResult = PushBack( reinterpret_cast<void*>(data), blockSize, kAliHLTDataTypeGlobalVertexer|kAliHLTDataOriginOut );  
+    fBenchmark.AddOutput(blockSize);
   }  
   
   
@@ -326,6 +318,7 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
     if( iResult==0 && data && data->fPrimNContributors >=3 ){
       AliESDVertex vESD( data->fPrimP, data->fPrimC, data->fPrimChi2, data->fPrimNContributors );
       iResult = PushBack( (TObject*) &vESD, kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut,0 );
+      fBenchmark.AddOutput(GetLastObjectSize());
     }
   }
 
@@ -333,24 +326,13 @@ int AliHLTGlobalVertexerComponent::DoEvent( const AliHLTComponentEventData& /*ev
   if( iResult==0 && event && data ){  
     FillESD( event, data ); 
     iResult = PushBack( event, kAliHLTDataTypeESDObject|kAliHLTDataOriginOut, 0);
+    fBenchmark.AddOutput(GetLastObjectSize());
   }
 
   delete[] buf;
 
-  timer.Stop();
-  fStatTimeR+=timer.RealTime();
-  fStatTimeC+=timer.CpuTime();
-  fStatNEvents++;
-
-  /*  
-  //if( fStatNEv%100==0 )
-  cout<<"SG: "<<GetComponentID()<<": "<<fStatNEvents<<" events, real time: total= "
-      <<fStatTimeR/fStatNEvents*1.e3<<" / create= "<<fStatTimeR1/fStatNEvents*1.e3
-      <<" / vprim= "<<fStatTimeR2/fStatNEvents*1.e3<<" / v0= "<<fStatTimeR3/fStatNEvents*1.e3
-      <<", CPU: total= "<<fStatTimeC/fStatNEvents*1.e3<<" / create= "<<fStatTimeC1/fStatNEvents*1.e3
-<<" / vprim= "<<fStatTimeC2/fStatNEvents*1.e3<<" / v0= "<<fStatTimeC2/fStatNEvents*1.e3
-      <<" ms"<<endl;
-  */
+  fBenchmark.Stop(0);
+  HLTInfo(fBenchmark.GetStatistics());
   return iResult;
 }
 
@@ -464,56 +446,114 @@ int AliHLTGlobalVertexerComponent::Reconfigure(const char* cdbEntry, const char*
 }
 
 
+struct AliHLTGlobalVertexerDeviation
+{
+  int fI; // track index
+  float fD; // deviation from vertex
+
+  bool operator<(const AliHLTGlobalVertexerDeviation &a) const { return fD<a.fD; }
+};
+
 
 void AliHLTGlobalVertexerComponent::FindPrimaryVertex()
 {
   //* Find event primary vertex
 
-  const AliKFParticle **vSelected = new const AliKFParticle*[fNTracks]; //* Selected particles for the vertex fit
-  Int_t *vIndex = new int [fNTracks];                    //* Indices of selected particles
-  Bool_t *vFlag = new bool [fNTracks];                    //* 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);
 
+  // select rough region (in sigmas) in which the vertex could be found, all tracks outside these limits are rejected
+  // from the primary vertex finding
   fPrimaryVtx.SetBeamConstraint( 0, 0, 0, 3., 3., 5.3 );
-  
+
+  const AliKFParticle **vSelected = new const AliKFParticle*[fNTracks]; //* Selected particles for the vertex fit
+  AliHLTGlobalVertexerDeviation *dev = new AliHLTGlobalVertexerDeviation[fNTracks];  
+
   Int_t nSelected = 0;
+  
   for( Int_t i = 0; i<fNTracks; 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;
+    dev[nSelected].fI = i; 
+    dev[nSelected].fD = chi; 
     vSelected[nSelected] = &(fTrackInfos[i].fParticle);
-    vIndex[nSelected] = i;
     nSelected++;  
+  }    
+
+  // fit
+
+  while( nSelected>2 ){ 
+
+    //* Primary vertex finder with rejection of outliers
+
+    for( Int_t i = 0; i<nSelected; i++){ 
+      vSelected[i] = &(fTrackInfos[dev[i].fI].fParticle);      
+    }
+
+    double xv = fPrimaryVtx.GetX();
+    double yv = fPrimaryVtx.GetY();
+    double zv = fPrimaryVtx.GetZ(); // values from previous iteration of calculations          
+    fPrimaryVtx.Initialize();
+    fPrimaryVtx.SetBeamConstraint( 0, 0, 0, 3., 3., 5.3 );
+    fPrimaryVtx.SetVtxGuess( xv, yv, zv );    
+
+    fPrimaryVtx.Construct( vSelected, nSelected, 0, -1, 1 ); // refilled for every iteration
+    
+    for( Int_t it=0; it<nSelected; it++ ){ 
+      const AliKFParticle &p = fTrackInfos[dev[it].fI].fParticle;
+      if( nSelected <= 20 ){
+       AliKFVertex tmp =  fPrimaryVtx - p; // exclude the current track from the sample and recalculate the vertex
+       dev[it].fD = p.GetDeviationFromVertex( tmp );
+      } else {
+       dev[it].fD = p.GetDeviationFromVertex( fPrimaryVtx );   
+      }
+    }
+    sort(dev,dev+nSelected); // sort tracks with increasing chi2 (used for rejection)
+       
+    int nRemove = (int) ( 0.3*nSelected );  //remove 30% of the tracks (done for performance, only if there are more than 20 tracks)  
+    if( nSelected - nRemove <=20 ) nRemove = 1;  // removal based on the chi2 of every track   
+    int firstRemove = nSelected - nRemove;
+    while( firstRemove<nSelected ){
+      if( dev[firstRemove].fD >= fConstrainedTrackDeviation ) break;
+      firstRemove++;
+    }
+    if( firstRemove>=nSelected ) break;
+    nSelected = firstRemove;
+  }
+
+  for( Int_t i = 0; i<fNTracks; i++){ 
+    fTrackInfos[i].fPrimUsedFlag = 0;
   }
-  fPrimaryVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, fConstrainedTrackDeviation );
 
+  if( nSelected < 3 ){  // no vertex for fewer than 3 contributors  
+    fPrimaryVtx.NDF() = -3;
+    fPrimaryVtx.Chi2() = 0;
+    nSelected = 0;
+  }
+  
   for( Int_t i = 0; i<nSelected; i++){ 
-    if( vFlag[i] ) fTrackInfos[vIndex[i]].fPrimUsedFlag = 1;
+    AliESDTrackInfo &info = fTrackInfos[dev[i].fI];
+    info.fPrimUsedFlag = 1;
+    info.fPrimDeviation = dev[i].fD;
   }
 
   for( Int_t i = 0; i<fNTracks; i++ ){
     AliESDTrackInfo &info = fTrackInfos[i];
+    if( info.fPrimUsedFlag ) continue;
     info.fPrimDeviation = info.fParticle.GetDeviationFromVertex( fPrimaryVtx );   
   }
 
-  //cout<<"SG: prim vtx event"<<fStatNEvents<<": n selected="<<nSelected<<", ncont="<<fPrimaryVtx.GetNContributors()<<endl;
-
-  if( fPrimaryVtx.GetNContributors()<3 ){
-    for( Int_t i = 0; i<fNTracks; i++)
-      fTrackInfos[i].fPrimUsedFlag = 0;
-  }
   delete[] vSelected;
-  delete[] vIndex;
-  delete[] vFlag;
+  delete[] dev;
 }
 
 
-void AliHLTGlobalVertexerComponent::FindV0s( vector<pair<int,int> > v0s  )
+
+void AliHLTGlobalVertexerComponent::FindV0s( vector<pair<int,int> > &v0s  )
 {
   //* V0 finder
 
@@ -543,7 +583,7 @@ void AliHLTGlobalVertexerComponent::FindV0s( vector<pair<int,int> > v0s  )
 
       if( (++statN)%100 ==0 ){ 
        if( timer.RealTime()>= fV0TimeLimit ){  run = 0; break; }
-       timer.Start();
+       timer.Continue();
       }
 
       //* check if the particles fit
@@ -623,6 +663,7 @@ void AliHLTGlobalVertexerComponent::FillESD( AliESDEvent *event, AliHLTGlobalVer
   if( data->fPrimNContributors >=3 ){
 
     AliESDVertex vESD( data->fPrimP, data->fPrimC, data->fPrimChi2, data->fPrimNContributors );
+    event->SetPrimaryVertexTPC( &vESD );
     event->SetPrimaryVertexTracks( &vESD );
 
     // relate tracks to the primary vertex