HLT/global/AliHLTPrimaryVertexFinderComponent.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE HLT Project        *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Timur Pocheptsov <Timur.Pocheptsov@cern.ch>           *
   8 //*                  for The ALICE HLT Project.                            *
   9 //*                                                                        *
  10 //* Permission to use, copy, modify and distribute this software and its   *
  11 //* documentation strictly for non-commercial purposes is hereby granted   *
  12 //* without fee, provided that the above copyright notice appears in all   *
  13 //* copies and that both the copyright notice and this permission notice   *
  14 //* appear in the supporting documentation. The authors make no claims     *
  15 //* about the suitability of this software for any purpose. It is          *
  16 //* provided "as is" without express or implied warranty.                  *
  17 //**************************************************************************
  18
  19 /// @file   AliHLTPrimaryVertexFinderComponent.cxx
  20 /// @author Timur Pocheptsov
  21 /// @date   2010-12-26
  22 /// @brief  Primary vertex finder component
  23 ///
  24
  25 #include <algorithm>
  26 #include <cerrno>
  27
  28 #include "TString.h"
  29 #include "TMath.h"
  30
  31 #include "AliHLTPrimaryVertexFinderComponent.h"
  32 #include "AliExternalTrackParam.h"
  33 #include "AliHLTDataTypes.h"
  34
  35 ClassImp(AliHLTPrimaryVertexFinderComponent)
  36
  37 const double AliHLTPrimaryVertexFinderComponent::fgDefaultDeviation = 4.;
  38
  39 //________________________________________________________________________
  40 AliHLTPrimaryVertexFinderComponent::AliHLTPrimaryVertexFinderComponent()
  41                              : fPrimaryOutput(),
  42                                fPrimaryVtx(),
  43                                fFitTracksToVertex(true),
  44                                fConstrainedTrackDeviation(fgDefaultDeviation)
  45 {
  46   //Default ctor.
  47 }
  48
  49 //________________________________________________________________________
  50 const char* AliHLTPrimaryVertexFinderComponent::GetComponentID()
  51 {
  52   //Component's "name".
  53   return "PrimaryVertexFinder";
  54 }
  55
  56 //________________________________________________________________________
  57 void AliHLTPrimaryVertexFinderComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
  58 {
  59   //Input to the primary vertex finder:
  60   //a)tracks from ESD object;
  61   //b)hlt tracks (ITS)
  62   //c)hlt tracks (TPC)
  63   list.clear();
  64
  65   list.push_back(kAliHLTDataTypeESDObject);
  66   list.push_back(kAliHLTDataTypeTrack | kAliHLTDataOriginITS);
  67   list.push_back(kAliHLTDataTypeTrack | kAliHLTDataOriginTPC);
  68 }
  69
  70 //________________________________________________________________________
  71 AliHLTComponentDataType AliHLTPrimaryVertexFinderComponent::GetOutputDataType()
  72 {
  73   //Data type of output.
  74   return kAliHLTMultipleDataType;
  75 }
  76
  77 //________________________________________________________________________
  78 int AliHLTPrimaryVertexFinderComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& list)
  79 {
  80   //Types for outputs from primary vertex finder (for V0 finder).
  81   list.clear();
  82
  83   //Indices of tracks, participating in a primary.
  84   list.push_back(kAliHLTDataTypePrimaryFinder | kAliHLTDataOriginOut);
  85   //KFVertex - primary vertex.
  86   list.push_back(kAliHLTDataTypeKFVertex | kAliHLTDataOriginOut);
  87
  88   return list.size();
  89 }
  90
  91 //________________________________________________________________________
  92 void AliHLTPrimaryVertexFinderComponent::GetOutputDataSize(unsigned long& constBase, double& inputMultiplier)
  93 {
  94   //These numbers are complete crap.
  95   constBase = 80000;
  96   inputMultiplier = 2.;
  97 }
  98
  99 //________________________________________________________________________
 100 AliHLTComponent* AliHLTPrimaryVertexFinderComponent::Spawn()
 101 {
 102   //Create primary vertex finder componet.
 103   return new AliHLTPrimaryVertexFinderComponent;
 104 }
 105
 106 //________________________________________________________________________
 107 int AliHLTPrimaryVertexFinderComponent::DoInit(int /*argc*/, const char** /*argv*/)
 108 {
 109   //Process options.
 110   //1. Default parameters.
 111   fFitTracksToVertex = true;
 112   fConstrainedTrackDeviation = fgDefaultDeviation;
 113
 114   //2. Parameters from OCDB.
 115   TString cdbPath("HLT/ConfigHLT/");
 116   cdbPath += GetComponentID();
 117
 118   //This part will be uncommented as soon as
 119   //OCDB object is added.
 120   /*
 121   int res = ConfigureFromCDBTObjString(cdbPath);
 122
 123   if (res < 0)
 124     return res;
 125
 126   //3. "Command line" parameters.
 127   if (argc)
 128     res = ConfigureFromArgumentString(argc, argv);
 129
 130   return res;*/
 131
 132   return 0;
 133 }
 134
 135 //________________________________________________________________________
 136 int AliHLTPrimaryVertexFinderComponent::ScanConfigurationArgument(int argc, const char** argv)
 137 {
 138   //Scan the name of option and its parameters from the list.
 139   //Return number of processed entries.
 140   //Possible arguments:
 141   //-fitTracksToVertex 1/0
 142   //-constrainedTrackDeviation value
 143   AliHLTUtility::CmdLineParser parser;
 144   parser.Add("-fitTrackToVertex", &fFitTracksToVertex);
 145   parser.Add("-constrainedTrackDeviation", &fConstrainedTrackDeviation);
 146
 147   const int nParsed = parser.Parse(argc, argv, 0);
 148   if (nParsed < 0) {
 149     HLTError(parser.GetError().Data());
 150     return -EPROTO;
 151   }
 152
 153   return nParsed;
 154 }
 155
 156 //________________________________________________________________________
 157 int AliHLTPrimaryVertexFinderComponent::DoDeinit()
 158 {
 159   //Reset parameters to default.
 160   fFitTracksToVertex = true;
 161   fConstrainedTrackDeviation = fgDefaultDeviation;
 162
 163   return 0;
 164 }
 165
 166 //________________________________________________________________________
 167 int AliHLTPrimaryVertexFinderComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/,
 168                                        AliHLTComponentTriggerData& /*trigData*/)
 169 {
 170   //Find primary vertex.
 171
 172   if (GetFirstInputBlock(kAliHLTDataTypeSOR) || GetFirstInputBlock(kAliHLTDataTypeEOR))
 173     return 0;
 174
 175
 176   //Clean all previous track infos and output block.
 177   fTrackInfos.clear();
 178   fPrimaryOutput.clear();
 179
 180   //Initialize KF package.
 181   AliKFParticle::SetField(GetBz());
 182
 183   //The logic, how vertex finder reads input tracks,
 184   //is taken from the original global vertexer.
 185   //First, try to read tracks from ESD event.
 186   //Both positive and negative PID are 211 ("default").
 187   //Other hypotesis are not possible here.
 188   ReadESDTracks(211, 211);
 189   //If no good esd tracks or no esd at all:
 190   if (!fTrackInfos.size())
 191     ReadHLTTracks(kAliHLTDataTypeTrack | kAliHLTDataOriginITS, 211, 211);
 192   //If no good its tracks:
 193   if (!fTrackInfos.size())
 194     ReadHLTTracks(kAliHLTDataTypeTrack | kAliHLTDataOriginTPC, 211, 211);
 195
 196   if (!fTrackInfos.size()) {
 197     HLTWarning("No input tracks found");
 198     return 0;
 199   }
 200
 201   FindPrimaryVertex();
 202
 203   return DoOutput();
 204 }
 205
 206 namespace
 207 {
 208
 209 struct VertexDeviation
 210 {
 211   int fI; //Index in fTrackInfos array.
 212   double fD; //Deviation from primary vertex.
 213
 214   bool operator < (const VertexDeviation& rhs) const
 215   {
 216     return fD < rhs.fD;
 217   }
 218 };
 219
 220 }
 221
 222 //________________________________________________________________________
 223 void AliHLTPrimaryVertexFinderComponent::FindPrimaryVertex()
 224 {
 225   //Find event's primary vertex.
 226
 227   ///////////////////////////////////////////////////////////////////////
 228   //Some changes must be done here to read the initial guess (?)
 229   //for primary vertex.
 230   //Select rough region (in sigmas) in which the vertex could be found,
 231   //all tracks outside these limits are rejected from the primary vertex finding.
 232   fPrimaryVtx.Initialize();
 233   fPrimaryVtx.SetBeamConstraint(0., 0., 0., 3., 3., 5.3);
 234   ////////////////////////////////////////////////////////////////////////
 235
 236   std::vector<const AliKFParticle*> vSelected(fTrackInfos.size());
 237   std::vector<VertexDeviation> devs(fTrackInfos.size());
 238
 239   int nSelected = 0;
 240   for (VectorSize_t i = 0; i < fTrackInfos.size(); ++i) {
 241     const AliKFParticle& p = fTrackInfos[i].fParticle;
 242     const double chi = p.GetDeviationFromVertex(fPrimaryVtx);
 243     if (chi > fConstrainedTrackDeviation)
 244       continue;
 245
 246     devs[nSelected].fI = i;
 247     devs[nSelected].fD = chi;
 248     vSelected[nSelected] = &fTrackInfos[i].fParticle;
 249     nSelected++;
 250   }
 251
 252   //Fit
 253   while (nSelected > 2) {
 254     //Primary vertex finder with rejection of outliers
 255     for (int i = 0; i < nSelected; ++i)
 256       vSelected[i] = &fTrackInfos[devs[i].fI].fParticle;
 257
 258     const double xv = fPrimaryVtx.GetX();
 259     const double yv = fPrimaryVtx.GetY();
 260     const double zv = fPrimaryVtx.GetZ(); //Values from the previous iteration.
 261
 262     fPrimaryVtx.Initialize();
 263     fPrimaryVtx.SetBeamConstraint(0, 0, 0, 3., 3., 5.3);
 264     fPrimaryVtx.SetVtxGuess(xv, yv, zv);
 265
 266     // refilled for every iteration
 267     //0: pointer to production vertex, -1. : mass, true : constrained.
 268     fPrimaryVtx.Construct(&vSelected[0], nSelected, 0, -1., true);
 269
 270     for (int it = 0; it < nSelected; ++it) {
 271       const AliKFParticle& p = fTrackInfos[devs[it].fI].fParticle;
 272       if (nSelected <= 20) {
 273         //Exclude the current track from the sample and recalculate the vertex
 274         AliKFVertex tmp(fPrimaryVtx - p);
 275         devs[it].fD = p.GetDeviationFromVertex(tmp);
 276       } else {
 277         devs[it].fD = p.GetDeviationFromVertex(fPrimaryVtx);
 278       }
 279     }
 280
 281     //Sort tracks with increasing chi2 (used for rejection)
 282     std::sort(&devs[0], &devs[0] + nSelected);
 283
 284     //Remove 30% of the tracks (done for performance, only if there are more than 20 tracks)
 285     int nRemove = int(0.3 * nSelected);
 286     if (nSelected - nRemove <= 20)
 287       nRemove = 1;// removal based on the chi2 of every track
 288
 289     int firstRemove = nSelected - nRemove;
 290     while (firstRemove < nSelected) {
 291       if (devs[firstRemove].fD >= fConstrainedTrackDeviation)
 292         break;
 293       firstRemove++;
 294     }
 295
 296     if (firstRemove >= nSelected)
 297       break;
 298
 299     nSelected = firstRemove;
 300   }
 301
 302   if (nSelected < 3) {//No vertex for less than 3 contributors.
 303     fPrimaryVtx.NDF() = -3;
 304     fPrimaryVtx.Chi2() = 0.;
 305     nSelected = 0;
 306   }
 307
 308   if (nSelected) {
 309     //Prepare output block.
 310     fPrimaryOutput.resize(sizeof(PrimaryFinderBlock) + sizeof(int) * (nSelected - 1));
 311     PrimaryFinderBlock* out = reinterpret_cast<PrimaryFinderBlock*>(&fPrimaryOutput[0]);
 312
 313     out->fFitTracksFlag = fFitTracksToVertex;
 314     out->fNPrimaryTracks = nSelected;
 315
 316     int minID = fTrackInfos[devs[0].fI].fID;
 317     int maxID = minID;
 318     for (int i = 0; i < nSelected; ++i) {
 319       const int id = fTrackInfos[devs[i].fI].fID;
 320       minID = TMath::Min(minID, id);
 321       maxID = TMath::Max(maxID, id);
 322       out->fPrimTrackIds[i] = id;
 323     }
 324
 325     out->fMinPrimID = minID;
 326     out->fMaxPrimID = maxID;
 327   }
 328 }
 329
 330 //________________________________________________________________________
 331 int AliHLTPrimaryVertexFinderComponent::DoOutput()
 332 {
 333   //Primary vertex finder output.
 334   if (!fPrimaryOutput.size()) {
 335     //Vertex not found.
 336     //Messages? return values?
 337     HLTWarning("No primary vertex was found");
 338     return 0;
 339   }
 340
 341   //1. indices of primary tracks;
 342   //int - type of PushBack's parameter.
 343   const int iResult = PushBack(&fPrimaryOutput[0], fPrimaryOutput.size(),
 344                                kAliHLTDataTypePrimaryFinder | kAliHLTDataOriginOut);
 345   if (iResult < 0)
 346     return iResult;
 347   //2. primary vertex (AliKFVertex).
 348   return PushBack(&fPrimaryVtx, kAliHLTDataTypeKFVertex | kAliHLTDataOriginOut, 0);
 349 }