Make SetChi2 method public
[u/mrichter/AliRoot.git] / HLT / global / AliHLTV0FinderComponent.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   AliHLTV0FinderComponent.cxx
20 /// @author Timur Pocheptsov
21 /// @date   2010-12-26
22 /// @brief  V0 finder component
23 ///
24
25 #include "TString.h"
26
27 #include "AliHLTDataTypes.h"
28 #include "AliHLTV0FinderComponent.h"
29
30 ClassImp(AliHLTV0FinderComponent)
31
32 const double AliHLTV0FinderComponent::fgDaughterPrimDeviation = 2.5;
33 const double AliHLTV0FinderComponent::fgPrimDeviation = 3.5;
34 const double AliHLTV0FinderComponent::fgChi = 3.5;
35 const double AliHLTV0FinderComponent::fgDecayLengthInSigmas = 3.;
36
37 //________________________________________________________________________
38 AliHLTV0FinderComponent::AliHLTV0FinderComponent()
39                   : fPrimaryVtx(),
40                     fPrimaryTracks(),
41                     fNPrimaryTracks(0),
42                     fMinPrimID(0),
43                     fMaxPrimID(0),
44                     fDaughterPrimDeviation(fgDaughterPrimDeviation),
45                     fPrimDeviation(fgPrimDeviation),
46                     fChi(fgChi),
47                     fDecayLengthInSigmas(fgDecayLengthInSigmas),
48                     fPosPID(211),
49                     fNegPID(211),
50                     fV0s(),
51                     fGammaFinder(false)
52 {
53   //Default ctor.
54 }
55
56 //________________________________________________________________________
57 const char* AliHLTV0FinderComponent::GetComponentID()
58 {
59   //Component's "name".
60   return "V0Finder";
61 }
62
63 //________________________________________________________________________
64 void AliHLTV0FinderComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
65 {
66   //Input to the primary vertex finder:
67   //a)tracks from ESD object;
68   //b)hlt tracks (ITS)
69   //c)hlt tracks (TPC)
70   //d)primary vertex (AliKFVertex)
71   //e)indices of primary tracks.
72   list.clear();
73   //Input tracks.
74   list.push_back(kAliHLTDataTypeESDObject);
75   list.push_back(kAliHLTDataTypeTrack | kAliHLTDataOriginITS);
76   list.push_back(kAliHLTDataTypeTrack | kAliHLTDataOriginTPC);
77   //Input from primary finder:
78   //Primary vertex.
79   list.push_back(kAliHLTDataTypeKFVertex | kAliHLTDataOriginOut);
80   //Primary tracks' indices.
81   list.push_back(kAliHLTDataTypePrimaryFinder | kAliHLTDataOriginOut);
82 }
83
84 //________________________________________________________________________
85 AliHLTComponentDataType AliHLTV0FinderComponent::GetOutputDataType()
86 {
87   //Data type of output.
88   return kAliHLTMultipleDataType;
89 }
90
91 //________________________________________________________________________
92 int AliHLTV0FinderComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
93 {
94   //Output type for V0 finder.
95   tgtList.clear();
96   //Indices of tracks, participating in V0s.
97   tgtList.push_back(kAliHLTDataTypeV0Finder | kAliHLTDataOriginOut);
98
99   return tgtList.size();
100 }
101
102 //________________________________________________________________________
103 void AliHLTV0FinderComponent::GetOutputDataSize(unsigned long& constBase, double& inputMultiplier)
104 {
105   //These numbers are complete crap.
106   constBase = 80000;
107   inputMultiplier = 2.;
108 }
109
110 //________________________________________________________________________
111 AliHLTComponent* AliHLTV0FinderComponent::Spawn()
112 {
113   //Create primary vertex finder componet.
114   return new AliHLTV0FinderComponent;
115 }
116
117 //________________________________________________________________________
118 int AliHLTV0FinderComponent::DoInit(int argc, const char** argv)
119 {
120   //1. Default parameters.
121   fDaughterPrimDeviation = fgDaughterPrimDeviation;
122   fPrimDeviation = fgPrimDeviation;
123   fChi = fgChi;
124   fDecayLengthInSigmas = fgDecayLengthInSigmas;
125   fPosPID = 211;
126   fNegPID = 211;
127   fGammaFinder = false;
128
129   //2. Parameters from OCDB.
130   TString cdbPath("HLT/ConfigHLT/");
131   cdbPath += GetComponentID();
132
133   int res = ConfigureFromCDBTObjString(cdbPath);
134   if (res < 0)
135     return res;
136
137   //3. "Command line" parameters.
138   if (argc)
139     res = ConfigureFromArgumentString(argc, argv);
140
141   fV0s.clear();
142   fV0s.push_back(0); //Number of v0s.
143
144   return res;
145 }
146
147 //________________________________________________________________________
148 int AliHLTV0FinderComponent::ScanConfigurationArgument(int argc, const char** argv)
149 {
150   //Scan one argument and its parameters from the list
151   //Return number of processed entries.
152   //Possible arguments:
153   //-daughterPrimDeviation num
154   //-primDeviation num
155   //-chi num
156   //-decayLengthInSigmas num
157   //-posPID int_num
158   //-negPid int_num
159   //-gammaFinder 0/1
160
161   AliHLTUtility::CmdLineParser parser;
162   parser.Add("-daughterPrimDeviation", &fDaughterPrimDeviation);
163   parser.Add("-primDeviation", &fPrimDeviation);
164   parser.Add("-chi", &fChi);
165   parser.Add("-decayLengthInSigmas", &fDecayLengthInSigmas);
166   parser.Add("-posPID", &fPosPID);
167   parser.Add("-negPID", &fNegPID);
168   parser.Add("-gammaFinder", &fGammaFinder);
169
170   const int nParsed = parser.Parse(argc, argv, 0);
171   if (nParsed < 0) {
172     HLTError(parser.GetError().Data());
173     return -EPROTO;
174   }
175
176   return nParsed;
177 }
178
179 //________________________________________________________________________
180 int AliHLTV0FinderComponent::DoDeinit()
181 {
182   //Reset parameters to default.
183   fDaughterPrimDeviation = fgDaughterPrimDeviation;
184   fPrimDeviation = fgPrimDeviation;
185   fChi = fgChi;
186   fDecayLengthInSigmas = fgDecayLengthInSigmas;
187   fPosPID = 211;
188   fNegPID = 211;
189   fGammaFinder = false;
190
191   return 0;
192 }
193
194 //________________________________________________________________________
195 int AliHLTV0FinderComponent::DoEvent(const AliHLTComponentEventData& /*evtData*/,
196                             AliHLTComponentTriggerData& /*trigData*/)
197 {
198   //Find primary vertex.
199   if (GetFirstInputBlock(kAliHLTDataTypeSOR) || GetFirstInputBlock(kAliHLTDataTypeEOR))
200     return 0;
201
202   //Clean all previous track infos.
203   fTrackInfos.clear();
204   fPrimaryTracks.clear();
205
206   //Initialize KF package.
207   AliKFParticle::SetField(GetBz());
208
209   // no output if there is no vertex information
210   if (!ReadPrimaryVertex())
211     return 0;
212
213   // start from clean array
214   // the output block contains at least the number of pairs even
215   // if this is zero
216   fV0s.clear();
217   fV0s.push_back(0); //Number of v0s.
218
219   if (ReadTracks()) {
220     FindV0s();
221   }
222
223   return DoOutput();
224 }
225
226 //________________________________________________________________________
227 bool AliHLTV0FinderComponent::ReadPrimaryVertex()
228 {
229   //Primary finder produces primary vertex (AliKFVertex) and
230   //indices for primary tracks.
231   //1. Try to extract primary vertex.
232   const TObject* obj = GetFirstInputObject(kAliHLTDataTypeKFVertex | kAliHLTDataOriginOut);
233   const AliKFVertex* kfVtx = dynamic_cast<const AliKFVertex*>(obj);
234
235   if (!kfVtx) {
236     HLTError("V0 finder requires KF vertex (primary vertex) as input");
237     return false;
238   }
239
240   //2. Try to read primary track indices.
241   const AliHLTComponentBlockData* p = GetFirstInputBlock(kAliHLTDataTypePrimaryFinder
242                                                           | kAliHLTDataOriginOut);
243   if (!p || !p->fSize || !p->fPtr) {
244     HLTError("Array of primary track indices expected");
245     return false;
246   }
247
248   //Data block from primary finder
249   const PrimaryFinderBlock* blk = static_cast<PrimaryFinderBlock*>(p->fPtr);
250   //Track ids must be positive integers.
251   if (blk->fMinPrimID < 0 || blk->fMaxPrimID < 0) {
252     HLTError("Got negative track ID from primary finder, internal HLT error");
253     return false;
254   }
255
256   //3. Got correct data, modify the component's state.
257   //KF vertex.
258   fPrimaryVtx = *kfVtx;
259   //Primary tracks.
260   fNPrimaryTracks = blk->fNPrimaryTracks;
261   fMinPrimID = blk->fMinPrimID;
262   fMaxPrimID = blk->fMaxPrimID;
263
264   fPrimaryTracks.assign(fMaxPrimID + 1, 0);
265   for (int i = 0; i < fNPrimaryTracks; ++i)
266     fPrimaryTracks[blk->fPrimTrackIds[i]] = 1;
267
268   return true;
269 }
270
271 //________________________________________________________________________
272 bool AliHLTV0FinderComponent::ReadTracks()
273 {
274   //The logic, how vertex finder reads input tracks,
275   //is taken from the original global vertexer.
276   //First, try to read tracks from ESD event.
277   ReadESDTracks(fPosPID, fNegPID);
278   if (fTrackInfos.size()) {
279     FindPrimaryDeviations();
280     return true;
281   }
282
283   //No good esd tracks, try:
284   ReadHLTTracks(kAliHLTDataTypeTrack | kAliHLTDataOriginITS, fPosPID, fNegPID);
285   if (fTrackInfos.size()) {
286     FindPrimaryDeviations();
287     return true;
288   }
289
290   //If no good its tracks, try:
291   ReadHLTTracks(kAliHLTDataTypeTrack | kAliHLTDataOriginTPC, fPosPID, fNegPID);
292   if (fTrackInfos.size()) {
293     FindPrimaryDeviations();
294     return true;
295   }
296
297   HLTError("No input tracks found for V0 finder");
298
299   return false;
300 }
301
302 //________________________________________________________________________
303 void AliHLTV0FinderComponent::FindPrimaryDeviations()
304 {
305   //Quite a tricky part.
306   for (VectorSize_t i = 0; i < fTrackInfos.size(); ++i) {
307     AliHLTTrackInfo& info = fTrackInfos[i];
308     if (IsPrimaryTrack(info.fID)) {
309       info.fPrimUsed = true;
310       //The way primary deviation is computed in primary finder:
311       if (fNPrimaryTracks <= 20) {
312         AliKFVertex tmp(fPrimaryVtx - info.fParticle);
313         info.fPrimDeviation = info.fParticle.GetDeviationFromVertex(tmp);
314       } else
315         info.fPrimDeviation = info.fParticle.GetDeviationFromVertex(fPrimaryVtx);
316     } else {
317       info.fPrimUsed = false;
318       info.fPrimDeviation = info.fParticle.GetDeviationFromVertex(fPrimaryVtx);
319     }
320   }
321 }
322
323 //________________________________________________________________________
324 bool AliHLTV0FinderComponent::IsPrimaryTrack(int id)const
325 {
326   if (id < fMinPrimID || id > fMaxPrimID)
327     return false;
328
329   return fPrimaryTracks[id];
330 }
331
332 //________________________________________________________________________
333 void AliHLTV0FinderComponent::FindV0s()
334 {
335   //Here's the core.
336   if (fPrimaryVtx.GetNContributors() < 3)
337     return;
338
339   for (int iTr = 0, ei = fTrackInfos.size(); iTr < ei; ++iTr) {
340     AliHLTTrackInfo& info = fTrackInfos[iTr];
341     if (info.fParticle.GetQ() > 0)
342       continue;
343     if (info.fPrimDeviation < fDaughterPrimDeviation)
344       continue;
345
346     for (int jTr = 0; jTr < ei; ++jTr) {
347       AliHLTTrackInfo& jnfo = fTrackInfos[jTr];
348       if (jnfo.fParticle.GetQ() < 0)
349         continue;
350       if (jnfo.fPrimDeviation < fDaughterPrimDeviation)
351         continue;
352
353       //Check if the particles fit
354       if (info.fParticle.GetDeviationFromParticle(jnfo.fParticle) > fChi)
355         continue;
356
357       //Construct V0 mother
358       AliKFParticle v0(info.fParticle, jnfo.fParticle /*, bGammaFinder*/);
359       //Check V0 Chi^2
360       if (v0.GetChi2() < 0. || v0.GetChi2() > fChi * fChi * v0.GetNDF())
361         continue;
362
363       //Subtruct daughters from primary vertex
364       AliKFVertex primVtxCopy(fPrimaryVtx);
365
366       if (info.fPrimUsed) {
367         if (primVtxCopy.GetNContributors() <= 2)
368           continue;
369         primVtxCopy -= info.fParticle;
370       }
371
372       if (jnfo.fPrimUsed) {
373         if (primVtxCopy.GetNContributors() <= 2)
374           continue;
375         primVtxCopy -= jnfo.fParticle;
376       }
377
378       //Check v0 Chi^2 deviation from primary vertex
379       if (v0.GetDeviationFromVertex(primVtxCopy) > fPrimDeviation)
380         continue;
381       //Add V0 to primary vertex to improve the primary vertex resolution
382       primVtxCopy += v0;
383       //Set production vertex for V0
384       v0.SetProductionVertex(primVtxCopy);
385       //Get V0 decay length with estimated error
386       double length = 0., sigmaLength = 0.;
387       if (v0.GetDecayLength(length, sigmaLength))
388         continue;
389       //Reject V0 if it decays too close[sigma] to the primary vertex
390       if (length  < fDecayLengthInSigmas * sigmaLength)
391         continue;
392       //Keep v0
393       fV0s.push_back(info.fID);
394       fV0s.push_back(jnfo.fID);
395
396       fV0s[0] += 1;
397     }
398   }
399 }
400
401 //________________________________________________________________________
402 int AliHLTV0FinderComponent::DoOutput()
403 {
404   //Save V0s' track pairs' indices.
405   if (!fV0s.size()) {
406     HLTError("internal data mismatch of output structure");
407     return 0;
408   }
409
410   HLTInfo("Number of v0 candidates: %d", fV0s[0]);
411
412   //Indices of primary tracks.
413   return PushBack(&fV0s[0], fV0s.size() * sizeof(int),
414                   kAliHLTDataTypeV0Finder | kAliHLTDataOriginOut);
415 }