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   ReadESDTracks();
 187   //If no good esd tracks or no esd at all:
 188   if (!fTrackInfos.size())
 189     ReadHLTTracks(kAliHLTDataTypeTrack | kAliHLTDataOriginITS);
 190   //If no good its tracks:
 191   if (!fTrackInfos.size())
 192     ReadHLTTracks(kAliHLTDataTypeTrack | kAliHLTDataOriginTPC);
 193
 194   if (!fTrackInfos.size()) {
 195     HLTWarning("No input tracks found");
 196     return 0;
 197   }
 198
 199   FindPrimaryVertex();
 200
 201   return DoOutput();
 202 }
 203
 204 namespace
 205 {
 206
 207 struct VertexDeviation
 208 {
 209   int fI; //Index in fTrackInfos array.
 210   double fD; //Deviation from primary vertex.
 211
 212   bool operator < (const VertexDeviation& rhs) const
 213   {
 214     return fD < rhs.fD;
 215   }
 216 };
 217
 218 }
 219
 220 //________________________________________________________________________
 221 void AliHLTPrimaryVertexFinderComponent::FindPrimaryVertex()
 222 {
 223   //Find event's primary vertex.
 224
 225   ///////////////////////////////////////////////////////////////////////
 226   //Some changes must be done here to read the initial guess (?)
 227   //for primary vertex.
 228   //Select rough region (in sigmas) in which the vertex could be found,
 229   //all tracks outside these limits are rejected from the primary vertex finding.
 230   fPrimaryVtx.Initialize();
 231   fPrimaryVtx.SetBeamConstraint(0., 0., 0., 3., 3., 5.3);
 232   ////////////////////////////////////////////////////////////////////////
 233
 234   std::vector<const AliKFParticle*> vSelected(fTrackInfos.size());
 235   std::vector<VertexDeviation> devs(fTrackInfos.size());
 236
 237   int nSelected = 0;
 238   for (VectorSize_t i = 0; i < fTrackInfos.size(); ++i) {
 239     const AliKFParticle& p = fTrackInfos[i].fParticle;
 240     const double chi = p.GetDeviationFromVertex(fPrimaryVtx);
 241     if (chi > fConstrainedTrackDeviation)
 242       continue;
 243
 244     devs[nSelected].fI = i;
 245     devs[nSelected].fD = chi;
 246     vSelected[nSelected] = &fTrackInfos[i].fParticle;
 247     nSelected++;
 248   }
 249
 250   //Fit
 251   while (nSelected > 2) {
 252     //Primary vertex finder with rejection of outliers
 253     for (int i = 0; i < nSelected; ++i)
 254       vSelected[i] = &fTrackInfos[devs[i].fI].fParticle;
 255
 256     const double xv = fPrimaryVtx.GetX();
 257     const double yv = fPrimaryVtx.GetY();
 258     const double zv = fPrimaryVtx.GetZ(); //Values from the previous iteration.
 259
 260     fPrimaryVtx.Initialize();
 261     fPrimaryVtx.SetBeamConstraint(0, 0, 0, 3., 3., 5.3);
 262     fPrimaryVtx.SetVtxGuess(xv, yv, zv);
 263
 264     // refilled for every iteration
 265     //0: pointer to production vertex, -1. : mass, true : constrained.
 266     fPrimaryVtx.Construct(&vSelected[0], nSelected, 0, -1., true);
 267
 268     for (int it = 0; it < nSelected; ++it) {
 269       const AliKFParticle& p = fTrackInfos[devs[it].fI].fParticle;
 270       if (nSelected <= 20) {
 271         //Exclude the current track from the sample and recalculate the vertex
 272         AliKFVertex tmp(fPrimaryVtx - p);
 273         devs[it].fD = p.GetDeviationFromVertex(tmp);
 274       } else {
 275         devs[it].fD = p.GetDeviationFromVertex(fPrimaryVtx);
 276       }
 277     }
 278
 279     //Sort tracks with increasing chi2 (used for rejection)
 280     std::sort(&devs[0], &devs[0] + nSelected);
 281
 282     //Remove 30% of the tracks (done for performance, only if there are more than 20 tracks)
 283     int nRemove = int(0.3 * nSelected);
 284     if (nSelected - nRemove <= 20)
 285       nRemove = 1;// removal based on the chi2 of every track
 286
 287     int firstRemove = nSelected - nRemove;
 288     while (firstRemove < nSelected) {
 289       if (devs[firstRemove].fD >= fConstrainedTrackDeviation)
 290         break;
 291       firstRemove++;
 292     }
 293
 294     if (firstRemove >= nSelected)
 295       break;
 296
 297     nSelected = firstRemove;
 298   }
 299
 300   if (nSelected < 3) {//No vertex for less than 3 contributors.
 301     fPrimaryVtx.NDF() = -3;
 302     fPrimaryVtx.Chi2() = 0.;
 303     nSelected = 0;
 304   }
 305
 306   if (nSelected) {
 307     //Prepare output block.
 308     fPrimaryOutput.resize(sizeof(PrimaryFinderBlock) + sizeof(int) * (nSelected - 1));
 309     PrimaryFinderBlock* out = reinterpret_cast<PrimaryFinderBlock*>(&fPrimaryOutput[0]);
 310
 311     out->fFitTracksFlag = fFitTracksToVertex;
 312     out->fNPrimaryTracks = nSelected;
 313
 314     int minID = fTrackInfos[devs[0].fI].fID;
 315     int maxID = minID;
 316     for (int i = 0; i < nSelected; ++i) {
 317       const int id = fTrackInfos[devs[i].fI].fID;
 318       minID = TMath::Min(minID, id);
 319       maxID = TMath::Max(maxID, id);
 320       out->fPrimTrackIds[i] = id;
 321     }
 322
 323     out->fMinPrimID = minID;
 324     out->fMaxPrimID = maxID;
 325   }
 326 }
 327
 328 //________________________________________________________________________
 329 int AliHLTPrimaryVertexFinderComponent::DoOutput()
 330 {
 331   //Primary vertex finder output.
 332   if (!fPrimaryOutput.size()) {
 333     //Vertex not found.
 334     //Messages? return values?
 335     HLTWarning("No primary vertex was found");
 336     return 0;
 337   }
 338
 339   //1. indices of primary tracks;
 340   //int - type of PushBack's parameter.
 341   const int iResult = PushBack(&fPrimaryOutput[0], fPrimaryOutput.size(),
 342                                kAliHLTDataTypePrimaryFinder | kAliHLTDataOriginOut);
 343   if (iResult < 0)
 344     return iResult;
 345   //2. primary vertex (AliKFVertex).
 346   return PushBack(&fPrimaryVtx, kAliHLTDataTypeKFVertex | kAliHLTDataOriginOut, 0);
 347 }