HLT/global/AliHLTVertexFinderBase.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   AliHLTVertexFinderBase.cxx
  20 /// @author Timur Pocheptsov
  21 /// @date   2010-12-26
  22 /// @brief  Base class for vertex finder components
  23 ///
  24
  25 #include <stdexcept>
  26 #include <cmath>
  27
  28 //AliHLTExternalTrackParam uses typedes from Rtypes.h
  29 //but does not include it.
  30 #include "Rtypes.h"
  31
  32 #include "AliHLTExternalTrackParam.h"
  33 #include "AliHLTGlobalBarrelTrack.h"
  34 #include "AliHLTVertexFinderBase.h"
  35 #include "AliExternalTrackParam.h"
  36 #include "AliESDVertex.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDtrack.h"
  39 #include "AliKFVertex.h"
  40 #include "AliESDv0.h"
  41
  42 ClassImp(AliHLTVertexFinderBase)
  43
  44 //_________________________________________________________________
  45 bool AliHLTVertexFinderBase::AliHLTTrackInfo::IsKFFinite()const
  46 {
  47   //Check KF particle.
  48
  49   //In C99, 'isfinite' is a macro.
  50   //But still, I add using directive, in case
  51   //it's a function in std.
  52   using namespace std;
  53
  54   for (int i = 0; i < 8; ++i)
  55     if (!isfinite(fParticle.GetParameter(i)))
  56       return false;
  57   for (int i = 0; i < 36; ++i)//Not sure, probably, 27 is enough.
  58     if (!isfinite(fParticle.GetCovariance(i)))
  59       return false;
  60
  61   return true;
  62 }
  63
  64 //_________________________________________________________________
  65 void AliHLTVertexFinderBase::ReadESDTracks(int posPID, int negPID)
  66 {
  67   //Try to read input tracks from AliESDEvent.
  68   //Use esd track's index as "track id" (fID).
  69   //IMPORTANT: fTrackInfos is _NOT_ cleared here,
  70   //must be done externally (if you need to).
  71
  72   const TObject* iter = GetFirstInputObject(kAliHLTDataTypeESDObject);
  73   for (; iter; iter = GetNextInputObject()) {
  74     //Cast away constness before dynamic_cast.
  75     AliESDEvent* esd = dynamic_cast<AliESDEvent*>((TObject*)iter);
  76     if (!esd)//From original code. Not sure, if this is possible at all.
  77       continue;
  78
  79     ReadESDTracks(esd, negPID, posPID);
  80     //Only one esd event is taken.
  81     break;
  82   }
  83 }
  84
  85 //_________________________________________________________________
  86 void AliHLTVertexFinderBase::ReadESDTracks(AliESDEvent* esd, int posPID, int negPID)
  87 {
  88   //Try to read input tracks from AliESDEvent.
  89   //Use esd track's index as "track id" (fID).
  90   //IMPORTANT: fTrackInfos is _NOT_ cleared here,
  91   //must be done externally (if you need to).
  92   esd->GetStdContent();
  93
  94   const int nTracks = esd->GetNumberOfTracks();
  95   if (nTracks)
  96     fTrackInfos.reserve(nTracks + fTrackInfos.size());
  97
  98   for (int i = 0; i < nTracks; ++i) {
  99     AliESDtrack* pTrack = esd->GetTrack(i);
 100     //This checks of track parameters are
 101     //from the original global vertexer.
 102     if (!pTrack)
 103       continue;
 104     if (pTrack->GetKinkIndex(0) > 0)
 105       continue;
 106     if (!(pTrack->GetStatus() & AliESDtrack::kTPCin))
 107       continue;
 108
 109     //i: track id, 211: pid, false: not used in primary, 0.: prim. deviation.
 110     AliHLTTrackInfo newTrackInfo(i, *pTrack, pTrack->Charge() > 0 ? posPID : negPID, false, 0.);
 111     if (!newTrackInfo.IsKFFinite())
 112       continue;
 113
 114     fTrackInfos.push_back(newTrackInfo);
 115   }
 116 }
 117
 118 //________________________________________________________________________
 119 void AliHLTVertexFinderBase::ReadHLTTracks(const AliHLTComponentDataType& blockType,
 120                                            int posPID, int negPID)
 121 {
 122   //Read HLT tracks as input.
 123   //IMPORTANT: fTrackInfos is _NOT_ cleared here,
 124   //must be done externally (if you need it).
 125
 126   const AliHLTComponentBlockData* block = GetFirstInputBlock(blockType);
 127   for (; block; block = GetNextInputBlock()) {
 128     const AliHLTTracksData* dataPtr = static_cast<AliHLTTracksData*>(block->fPtr);
 129     const AliHLTExternalTrackParam* hltTrk = dataPtr->fTracklets;
 130
 131     const int nTracks = dataPtr->fCount;
 132     if (nTracks)
 133       fTrackInfos.reserve(fTrackInfos.size() + nTracks);
 134
 135     for (int i = 0; i < nTracks; ++i) {
 136       //Ugly conversion from one track format to another
 137       //and to AliKFParitcle then, to be changed later.
 138       AliHLTGlobalBarrelTrack tmpTrk(*hltTrk);
 139
 140       AliHLTTrackInfo newTrackInfo(hltTrk->fTrackID, tmpTrk,
 141                                    tmpTrk.Charge() > 0 ? posPID : negPID, false, 0.);
 142       if (!newTrackInfo.IsKFFinite())
 143         continue;
 144
 145       fTrackInfos.push_back(newTrackInfo);
 146
 147       const unsigned dSize = sizeof(AliHLTExternalTrackParam)
 148                             + hltTrk->fNPoints * sizeof(unsigned);
 149       hltTrk = (AliHLTExternalTrackParam*)((char *)hltTrk + dSize);
 150     }
 151   }
 152 }
 153
 154 //________________________________________________________________________
 155 void AliHLTVertexFinderBase::FillESD(AliESDEvent* esd, AliKFVertex* vtx,
 156                                      const void* primData, const void* v0Data)
 157 {
 158   //Put the output of a primary finder and v0 finder to the esd event.
 159   //Code was taken from the original global vertexer.
 160   //Code is an absolute mess.
 161
 162   typedef std::map<int, int> Map_t;
 163   typedef Map_t::const_iterator MapIter_t;
 164
 165   //Map esdId -> esdTrackIndex.
 166   Map_t mapId;
 167
 168   double params[3];
 169   for (int i = 0; i < 3; ++i)
 170     params[i] = vtx->Parameters()[i];
 171   double cov[6];
 172   for (int i = 0; i < 6; ++i)
 173     cov[i] = vtx->CovarianceMatrix()[i];
 174
 175   AliESDVertex vESD(params, cov, vtx->GetChi2(), vtx->GetNContributors());
 176   esd->SetPrimaryVertexTPC(&vESD);
 177   esd->SetPrimaryVertexTracks(&vESD);
 178
 179   //Relate tracks to the primary vertex
 180   const PrimaryFinderBlock* prim = static_cast<const PrimaryFinderBlock*>(primData);
 181   const int nESDTracks = esd->GetNumberOfTracks();
 182   std::vector<bool> constrainedToVtx(nESDTracks);
 183
 184   if (prim->fFitTracksFlag) {
 185     for (int i = 0; i < nESDTracks; ++i) {
 186       if (!esd->GetTrack(i))
 187         continue;
 188
 189       mapId[esd->GetTrack(i)->GetID()] = i;
 190     }
 191
 192     for (int i = 0; i < prim->fNPrimaryTracks; ++i) {
 193       MapIter_t it = mapId.find(prim->fPrimTrackIds[i]);
 194       if (it == mapId.end())
 195         continue;
 196       const int itr = it->second;
 197       //100. is an argument for parameter maxd in AliESDtrack - cut on impact parameter.
 198       esd->GetTrack(itr)->RelateToVertex(&vESD, esd->GetMagneticField(), 100.);
 199       constrainedToVtx[itr] = true;
 200     }
 201   }
 202
 203   //Add v0s.
 204   const int* v0s = static_cast<const int*>(v0Data);
 205   const int nV0s = v0s[0];
 206   ++v0s;
 207   for (int i = 0; i < nV0s; ++i) {
 208     MapIter_t it = mapId.find(v0s[2 * i]);
 209     if (it==mapId.end())
 210       continue;
 211     const int iTr = it->second;
 212
 213     it = mapId.find(v0s[2 * i + 1]);
 214     if (it == mapId.end())
 215       continue;
 216     const int jTr = it->second;
 217
 218     AliESDv0 v0(*esd->GetTrack(iTr), iTr, *esd->GetTrack(jTr), jTr);
 219     esd->AddV0(&v0);
 220     // relate the tracks to the vertex
 221     if (prim->fFitTracksFlag) {
 222       if (constrainedToVtx[iTr] || constrainedToVtx[jTr])
 223         continue;
 224
 225       double pos[3];
 226       double sigma[3] = {.1, .1, .1};
 227       v0.XvYvZv(pos);
 228       AliESDVertex v0ESD(pos, sigma);
 229       esd->GetTrack(iTr)->RelateToVertex(&v0ESD, esd->GetMagneticField(), 100.);
 230       esd->GetTrack(jTr)->RelateToVertex(&v0ESD, esd->GetMagneticField(), 100.);
 231       constrainedToVtx[iTr] = true;
 232       constrainedToVtx[jTr] = true;
 233     }
 234   }
 235 }
 236
 237 namespace AliHLTUtility
 238 {
 239
 240 //_______________________________________________________________________
 241 Parameter::Parameter()
 242              : fWasSet(false),
 243                fConstraint(none),
 244                fBool(0),
 245                fInt(0),
 246                fDouble(0),
 247                fCompound(0),
 248                fConstraintChecker(0)
 249 {
 250   //Default ctor.
 251 }
 252
 253 //_______________________________________________________________________
 254 Parameter::Parameter(bool* b)
 255              : fWasSet(false),
 256                fConstraint(none),
 257                fBool(b),
 258                fInt(0),
 259                fDouble(0),
 260                fCompound(0),
 261                fConstraintChecker(0)
 262 {
 263   //Parameter of type bool.
 264 }
 265
 266 //_______________________________________________________________________
 267 Parameter::Parameter(int* i, Constraint c)
 268              : fWasSet(false),
 269                fConstraint(c),
 270                fBool(0),
 271                fInt(i),
 272                fDouble(0),
 273                fCompound(0),
 274                fConstraintChecker(0)
 275 {
 276   //Parameter of type int.
 277   SetConstraintChecker();
 278 }
 279
 280 //_______________________________________________________________________
 281 Parameter::Parameter(double* d, Constraint c)
 282              : fWasSet(false),
 283                fConstraint(c),
 284                fBool(0),
 285                fInt(0),
 286                fDouble(d),
 287                fCompound(0),
 288                fConstraintChecker(0)
 289 {
 290   //Parameter of type double.
 291   SetConstraintChecker();
 292 }
 293
 294 //_______________________________________________________________________
 295 Parameter::Parameter(CompoundType *ct)
 296              : fWasSet(false),
 297                fConstraint(none),
 298                fBool(0),
 299                fInt(0),
 300                fDouble(0),
 301                fCompound(ct),
 302                fConstraintChecker(0)
 303 {
 304   //Parameter of more complex user-defined type.
 305   //All checks and conversions must be implemented
 306   //by user of CompoundType.
 307 }
 308
 309 //_______________________________________________________________________
 310 unsigned Parameter::SetParameter(unsigned argc, const char** argv, unsigned currPos)
 311 {
 312   //Set parameter from command line tokens.
 313
 314   //It's up to compound parameter to parse.
 315   if (fCompound)
 316     return fCompound->SetParameter(argc, argv, currPos);
 317
 318   //Now, int, bool or double must be set from a string.
 319   //Primitive checks are done here.
 320   if (currPos == argc)
 321     throw std::runtime_error("value expected");
 322   if (fWasSet)
 323     throw std::runtime_error("parameter was set already");
 324
 325   const TString val(argv[currPos]);
 326   if (!val.Length())
 327     throw std::runtime_error("value expected");
 328
 329   if (fBool) {
 330     if (val.Length() != 1 || (val[0] != '0' && val[0] != '1'))
 331       throw std::runtime_error("expected 0 or 1 for bool parameter");
 332     *fBool = val[0] - '0';
 333   } else if (fInt)
 334     *fInt = val.Atoi();
 335   else
 336     *fDouble = val.Atof();
 337
 338   if (fConstraintChecker)
 339     fConstraintChecker(*this);
 340
 341   fWasSet = true;
 342
 343   //For double, int and bool - only one element of argv is processed.
 344   return 1;
 345 }
 346
 347 //_______________________________________________________________________
 348 void Parameter::SetConstraintChecker()
 349 {
 350   //Set function pointer.
 351   if (fConstraint == positive)
 352     fConstraintChecker = CheckPositive;
 353   else if (fConstraint == nonNegative)
 354     fConstraintChecker = CheckNonNegative;
 355   //No constraints for bool or CompoundType.
 356 }
 357
 358 //Aux. functions to check constraints.
 359 //_______________________________________________________________________
 360 void Parameter::CheckPositive(const Parameter& param)
 361 {
 362   //val > 0
 363   if (param.fInt) {
 364     if (*param.fInt <= 0)
 365       throw std::runtime_error("integer parameter must be positive");
 366   } else if (param.fDouble) {
 367      if (*param.fDouble <= 0.)//Hmmmm.
 368       throw std::runtime_error("double parameter must be positive");
 369   }
 370   //For bool or CompoundType there is no check.
 371 }
 372
 373 //_______________________________________________________________________
 374 void Parameter::CheckNonNegative(const Parameter& param)
 375 {
 376   //val >= 0
 377   if (param.fInt) {
 378     if (*param.fInt < 0)
 379       throw std::runtime_error("integer parameter must be non-negative");
 380   } else if (param.fDouble) {
 381     if (*param.fDouble < 0.)
 382       throw std::runtime_error("double parameter must be non-negative");
 383   }
 384   //For bool or CompoundType there is no check.
 385 }
 386
 387
 388 //Command line arguments parser.
 389
 390 //_______________________________________________________________________
 391 void CmdLineParser::Add(const TString& cmd, bool* b)
 392 {
 393   //Command with bool parameter.
 394   if (fParameters.find(cmd) == fParameters.end())
 395     fParameters[cmd] = Parameter(b);
 396 }
 397
 398 //_______________________________________________________________________
 399 void CmdLineParser::Add(const TString& cmd, int* i, Parameter::Constraint c)
 400 {
 401   //Command with int parameter, with constraint.
 402   if (fParameters.find(cmd) == fParameters.end())
 403     fParameters[cmd] = Parameter(i, c);
 404 }
 405
 406 //_______________________________________________________________________
 407 void CmdLineParser::Add(const TString& cmd, double* d, Parameter::Constraint c)
 408 {
 409   //Coomand with a parameter of type double, possibly with a constraint.
 410   if (fParameters.find(cmd) == fParameters.end())
 411     fParameters[cmd] = Parameter(d, c);
 412 }
 413
 414 //_______________________________________________________________________
 415 void CmdLineParser::Add(const TString& cmd, CompoundType *ct)
 416 {
 417   //Command with a parameter of compound type.
 418   if (fParameters.find(cmd) == fParameters.end())
 419     fParameters[cmd] = Parameter(ct);
 420 }
 421
 422 //_______________________________________________________________________
 423 void CmdLineParser::IgnoreCommand(const TString& cmd, unsigned nParams)
 424 {
 425   //Command and number of parameters to ignore.
 426   //I do not check, if it's in map already.
 427   fIgnored[cmd] = nParams;
 428 }
 429
 430 //_______________________________________________________________________
 431 int CmdLineParser::Parse(unsigned argc, const char** argv, unsigned currPos)
 432 {
 433   //Parse. Returns number of processed arguments.
 434   //Negative number - error (as it's in HLT)
 435   fError.Clear();
 436
 437   unsigned nProcessed = 0;
 438   while (currPos < argc) {
 439     //Command.
 440     const TString command(argv[currPos]);
 441     //Check, if this command must be ignored.
 442     SkipMapIter_t skipCmd = fIgnored.find(command);
 443     if (skipCmd != fIgnored.end()) {
 444       //Command and its parameters must be skipped
 445       nProcessed += 1 + skipCmd->second;
 446       currPos += 1 + skipCmd->second;
 447       continue;
 448     }
 449
 450     //Corresponding parameter.
 451     MapIter_t it = fParameters.find(command);
 452     if (it == fParameters.end()) {
 453       fError = "Unknown command: " + command;
 454       return -1;
 455     }
 456
 457     ++currPos;
 458     ++nProcessed;
 459
 460     try {
 461       const unsigned n = it->second.SetParameter(argc, argv, currPos);
 462       nProcessed += n;
 463       currPos += n;
 464     } catch (const std::runtime_error& e) {
 465       fError = "Command " + command + " error: " + e.what();
 466       return -1;
 467     }
 468   }
 469
 470   return nProcessed;
 471 }
 472
 473 }