initializing variables that were declare without a value
[u/mrichter/AliRoot.git] / 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 }