Make SetChi2 method public
[u/mrichter/AliRoot.git] / 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 }