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   if (v0Data) {
 205   const int* v0s = static_cast<const int*>(v0Data);
 206   const int nV0s = v0s[0];
 207   ++v0s;
 208   for (int i = 0; i < nV0s; ++i) {
 209     MapIter_t it = mapId.find(v0s[2 * i]);
 210     if (it==mapId.end())
 211       continue;
 212     const int iTr = it->second;
 213
 214     it = mapId.find(v0s[2 * i + 1]);
 215     if (it == mapId.end())
 216       continue;
 217     const int jTr = it->second;
 218
 219     AliESDv0 v0(*esd->GetTrack(iTr), iTr, *esd->GetTrack(jTr), jTr);
 220     esd->AddV0(&v0);
 221     // relate the tracks to the vertex
 222     if (prim->fFitTracksFlag) {
 223       if (constrainedToVtx[iTr] || constrainedToVtx[jTr])
 224         continue;
 225
 226       double pos[3];
 227       double sigma[3] = {.1, .1, .1};
 228       v0.XvYvZv(pos);
 229       AliESDVertex v0ESD(pos, sigma);
 230       esd->GetTrack(iTr)->RelateToVertex(&v0ESD, esd->GetMagneticField(), 100.);
 231       esd->GetTrack(jTr)->RelateToVertex(&v0ESD, esd->GetMagneticField(), 100.);
 232       constrainedToVtx[iTr] = true;
 233       constrainedToVtx[jTr] = true;
 234     }
 235   }
 236   }
 237 }
 238
 239 namespace AliHLTUtility
 240 {
 241
 242 //_______________________________________________________________________
 243 Parameter::Parameter()
 244              : fWasSet(false),
 245                fConstraint(none),
 246                fBool(0),
 247                fInt(0),
 248                fDouble(0),
 249                fCompound(0),
 250                fConstraintChecker(0)
 251 {
 252   //Default ctor.
 253 }
 254
 255 //_______________________________________________________________________
 256 Parameter::Parameter(bool* b)
 257              : fWasSet(false),
 258                fConstraint(none),
 259                fBool(b),
 260                fInt(0),
 261                fDouble(0),
 262                fCompound(0),
 263                fConstraintChecker(0)
 264 {
 265   //Parameter of type bool.
 266 }
 267
 268 //_______________________________________________________________________
 269 Parameter::Parameter(int* i, Constraint c)
 270              : fWasSet(false),
 271                fConstraint(c),
 272                fBool(0),
 273                fInt(i),
 274                fDouble(0),
 275                fCompound(0),
 276                fConstraintChecker(0)
 277 {
 278   //Parameter of type int.
 279   SetConstraintChecker();
 280 }
 281
 282 //_______________________________________________________________________
 283 Parameter::Parameter(double* d, Constraint c)
 284              : fWasSet(false),
 285                fConstraint(c),
 286                fBool(0),
 287                fInt(0),
 288                fDouble(d),
 289                fCompound(0),
 290                fConstraintChecker(0)
 291 {
 292   //Parameter of type double.
 293   SetConstraintChecker();
 294 }
 295
 296 //_______________________________________________________________________
 297 Parameter::Parameter(CompoundType *ct)
 298              : fWasSet(false),
 299                fConstraint(none),
 300                fBool(0),
 301                fInt(0),
 302                fDouble(0),
 303                fCompound(ct),
 304                fConstraintChecker(0)
 305 {
 306   //Parameter of more complex user-defined type.
 307   //All checks and conversions must be implemented
 308   //by user of CompoundType.
 309 }
 310
 311 //_______________________________________________________________________
 312 unsigned Parameter::SetParameter(unsigned argc, const char** argv, unsigned currPos)
 313 {
 314   //Set parameter from command line tokens.
 315
 316   //It's up to compound parameter to parse.
 317   if (fCompound)
 318     return fCompound->SetParameter(argc, argv, currPos);
 319
 320   //Now, int, bool or double must be set from a string.
 321   //Primitive checks are done here.
 322   if (currPos == argc)
 323     throw std::runtime_error("value expected");
 324   if (fWasSet)
 325     throw std::runtime_error("parameter was set already");
 326
 327   const TString val(argv[currPos]);
 328   if (!val.Length())
 329     throw std::runtime_error("value expected");
 330
 331   if (fBool) {
 332     if (val.Length() != 1 || (val[0] != '0' && val[0] != '1'))
 333       throw std::runtime_error("expected 0 or 1 for bool parameter");
 334     *fBool = val[0] - '0';
 335   } else if (fInt)
 336     *fInt = val.Atoi();
 337   else
 338     *fDouble = val.Atof();
 339
 340   if (fConstraintChecker)
 341     fConstraintChecker(*this);
 342
 343   fWasSet = true;
 344
 345   //For double, int and bool - only one element of argv is processed.
 346   return 1;
 347 }
 348
 349 //_______________________________________________________________________
 350 void Parameter::SetConstraintChecker()
 351 {
 352   //Set function pointer.
 353   if (fConstraint == positive)
 354     fConstraintChecker = CheckPositive;
 355   else if (fConstraint == nonNegative)
 356     fConstraintChecker = CheckNonNegative;
 357   //No constraints for bool or CompoundType.
 358 }
 359
 360 //Aux. functions to check constraints.
 361 //_______________________________________________________________________
 362 void Parameter::CheckPositive(const Parameter& param)
 363 {
 364   //val > 0
 365   if (param.fInt) {
 366     if (*param.fInt <= 0)
 367       throw std::runtime_error("integer parameter must be positive");
 368   } else if (param.fDouble) {
 369      if (*param.fDouble <= 0.)//Hmmmm.
 370       throw std::runtime_error("double parameter must be positive");
 371   }
 372   //For bool or CompoundType there is no check.
 373 }
 374
 375 //_______________________________________________________________________
 376 void Parameter::CheckNonNegative(const Parameter& param)
 377 {
 378   //val >= 0
 379   if (param.fInt) {
 380     if (*param.fInt < 0)
 381       throw std::runtime_error("integer parameter must be non-negative");
 382   } else if (param.fDouble) {
 383     if (*param.fDouble < 0.)
 384       throw std::runtime_error("double parameter must be non-negative");
 385   }
 386   //For bool or CompoundType there is no check.
 387 }
 388
 389
 390 //Command line arguments parser.
 391
 392 //_______________________________________________________________________
 393 void CmdLineParser::Add(const TString& cmd, bool* b)
 394 {
 395   //Command with bool parameter.
 396   if (fParameters.find(cmd) == fParameters.end())
 397     fParameters[cmd] = Parameter(b);
 398 }
 399
 400 //_______________________________________________________________________
 401 void CmdLineParser::Add(const TString& cmd, int* i, Parameter::Constraint c)
 402 {
 403   //Command with int parameter, with constraint.
 404   if (fParameters.find(cmd) == fParameters.end())
 405     fParameters[cmd] = Parameter(i, c);
 406 }
 407
 408 //_______________________________________________________________________
 409 void CmdLineParser::Add(const TString& cmd, double* d, Parameter::Constraint c)
 410 {
 411   //Coomand with a parameter of type double, possibly with a constraint.
 412   if (fParameters.find(cmd) == fParameters.end())
 413     fParameters[cmd] = Parameter(d, c);
 414 }
 415
 416 //_______________________________________________________________________
 417 void CmdLineParser::Add(const TString& cmd, CompoundType *ct)
 418 {
 419   //Command with a parameter of compound type.
 420   if (fParameters.find(cmd) == fParameters.end())
 421     fParameters[cmd] = Parameter(ct);
 422 }
 423
 424 //_______________________________________________________________________
 425 void CmdLineParser::IgnoreCommand(const TString& cmd, unsigned nParams)
 426 {
 427   //Command and number of parameters to ignore.
 428   //I do not check, if it's in map already.
 429   fIgnored[cmd] = nParams;
 430 }
 431
 432 //_______________________________________________________________________
 433 int CmdLineParser::Parse(unsigned argc, const char** argv, unsigned currPos)
 434 {
 435   //Parse. Returns number of processed arguments.
 436   //Negative number - error (as it's in HLT)
 437   fError.Clear();
 438
 439   unsigned nProcessed = 0;
 440   while (currPos < argc) {
 441     //Command.
 442     const TString command(argv[currPos]);
 443     //Check, if this command must be ignored.
 444     SkipMapIter_t skipCmd = fIgnored.find(command);
 445     if (skipCmd != fIgnored.end()) {
 446       //Command and its parameters must be skipped
 447       nProcessed += 1 + skipCmd->second;
 448       currPos += 1 + skipCmd->second;
 449       continue;
 450     }
 451
 452     //Corresponding parameter.
 453     MapIter_t it = fParameters.find(command);
 454     if (it == fParameters.end()) {
 455       fError = "Unknown command: " + command;
 456       return -1;
 457     }
 458
 459     ++currPos;
 460     ++nProcessed;
 461
 462     try {
 463       const unsigned n = it->second.SetParameter(argc, argv, currPos);
 464       nProcessed += n;
 465       currPos += n;
 466     } catch (const std::runtime_error& e) {
 467       fError = "Command " + command + " error: " + e.what();
 468       return -1;
 469     }
 470   }
 471
 472   return nProcessed;
 473 }
 474
 475 }