AliPhysics  a6017e1 (a6017e1)
AliAnalysisTaskSEImproveITS.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 
16 #include <TObjArray.h>
17 #include <TClonesArray.h>
18 #include <TGraph.h>
19 #include <TFile.h>
20 #include <TList.h>
21 #include <TNtuple.h>
22 
23 #include "AliVertex.h"
24 #include "AliVVertex.h"
25 #include "AliESDVertex.h"
26 #include "AliESDEvent.h"
27 #include "AliVertexerTracks.h"
28 #include "AliAODEvent.h"
29 #include "AliAODTrack.h"
30 #include "AliMCEvent.h"
31 #include "AliAODMCParticle.h"
32 #include "AliExternalTrackParam.h"
35 #include "AliAODRecoCascadeHF.h"
36 #include "AliAnalysisVertexingHF.h"
37 #include "AliNeutralTrackParam.h"
39 
40 //
41 // Implementation of the "hybrid-approach" for ITS upgrade studies.
42 // The tastk smears the track parameters according to estimations
43 // from single-track upgrade studies. Afterwards it recalculates
44 // the parameters of the reconstructed decays.
45 //
46 // WARNING: This will affect all tasks in a train after this one
47 // (which is typically desired, though).
48 //
49 
52  fD0ZResPCur (0),
53  fD0ZResKCur (0),
54  fD0ZResPiCur (0),
55  fD0ZResECur (0),
56  fD0RPResPCur (0),
57  fD0RPResKCur (0),
58  fD0RPResPiCur(0),
59  fD0RPResECur(0),
60  fD0RPSigmaPullRatioP(0),
61  fD0RPSigmaPullRatioK(0),
62  fD0RPSigmaPullRatioPi(0),
63  fD0RPSigmaPullRatioE(0),
64  fPt1ResPCur (0),
65  fPt1ResKCur (0),
66  fPt1ResPiCur (0),
67  fPt1ResECur (0),
68  fD0ZResPUpg (0),
69  fD0ZResKUpg (0),
70  fD0ZResPiUpg (0),
71  fD0ZResEUpg (0),
72  fD0RPResPUpg (0),
73  fD0RPResKUpg (0),
74  fD0RPResPiUpg(0),
75  fD0RPResEUpg(0),
76  fPt1ResPUpg (0),
77  fPt1ResKUpg (0),
78  fPt1ResPiUpg (0),
79  fPt1ResEUpg (0),
80  fD0ZResPCurSA (0),
81  fD0ZResKCurSA (0),
82  fD0ZResPiCurSA (0),
83  fD0ZResECurSA (0),
84  fD0RPResPCurSA (0),
85  fD0RPResKCurSA (0),
86  fD0RPResPiCurSA(0),
87  fD0RPResECurSA(0),
88  fPt1ResPCurSA (0),
89  fPt1ResKCurSA (0),
90  fPt1ResPiCurSA (0),
91  fPt1ResECurSA (0),
92  fD0ZResPUpgSA (0),
93  fD0ZResKUpgSA (0),
94  fD0ZResPiUpgSA (0),
95  fD0ZResEUpgSA (0),
96  fD0RPResPUpgSA (0),
97  fD0RPResKUpgSA (0),
98  fD0RPResPiUpgSA(0),
99  fD0RPResEUpgSA(0),
100  fPt1ResPUpgSA (0),
101  fPt1ResKUpgSA (0),
102  fPt1ResPiUpgSA (0),
103  fPt1ResEUpgSA (0),
104  fRunInVertexing(kFALSE),
105  fImproveTracks(kTRUE),
106  fUpdateSecVertCovMat(kTRUE),
107  fUpdateSTCovMatrix(kTRUE),
108  fUpdatePulls(kTRUE),
109  fMimicData(kFALSE),
110  fIsAOD (kTRUE),
111  fMCs (0),
112  fDebugOutput (0),
113  fDebugNtuple (0),
114  fDebugVars (0),
115  fNDebug (0)
116 {
117  //
118  // Default constructor.
119  for(Int_t jh=0; jh<2; jh++){
120  for(Int_t ih=0; ih<4; ih++){
121  fD0RPMeanPCur[jh][ih]=0x0;
122  fD0RPMeanKCur[jh][ih]=0x0;
123  fD0RPMeanPiCur[jh][ih]=0x0;
124  fD0RPMeanECur[jh][ih]=0x0;
125  fD0RPMeanPUpg[jh][ih]=0x0;
126  fD0RPMeanKUpg[jh][ih]=0x0;
127  fD0RPMeanPiUpg[jh][ih]=0x0;
128  fD0RPMeanEUpg[jh][ih]=0x0;
129  }
130  }
131  //
132 }
133 
135  const char *period,
136  const char *systematic,
137  Bool_t isRunInVertexing,
138  Int_t ndebug)
139  :AliAnalysisTaskSE(name),
140  fD0ZResPCur (0),
141  fD0ZResKCur (0),
142  fD0ZResPiCur (0),
143  fD0ZResECur (0),
144  fD0RPResPCur (0),
145  fD0RPResKCur (0),
146  fD0RPResPiCur(0),
147  fD0RPResECur(0),
152  fPt1ResPCur (0),
153  fPt1ResKCur (0),
154  fPt1ResPiCur (0),
155  fPt1ResECur (0),
156  fD0ZResPUpg (0),
157  fD0ZResKUpg (0),
158  fD0ZResPiUpg (0),
159  fD0ZResEUpg (0),
160  fD0RPResPUpg (0),
161  fD0RPResKUpg (0),
162  fD0RPResPiUpg(0),
163  fD0RPResEUpg(0),
164  fPt1ResPUpg (0),
165  fPt1ResKUpg (0),
166  fPt1ResPiUpg (0),
167  fPt1ResEUpg (0),
168  fD0ZResPCurSA (0),
169  fD0ZResKCurSA (0),
170  fD0ZResPiCurSA (0),
171  fD0ZResECurSA (0),
172  fD0RPResPCurSA (0),
173  fD0RPResKCurSA (0),
174  fD0RPResPiCurSA(0),
175  fD0RPResECurSA(0),
176  fPt1ResPCurSA (0),
177  fPt1ResKCurSA (0),
178  fPt1ResPiCurSA (0),
179  fPt1ResECurSA (0),
180  fD0ZResPUpgSA (0),
181  fD0ZResKUpgSA (0),
182  fD0ZResPiUpgSA (0),
183  fD0ZResEUpgSA (0),
184  fD0RPResPUpgSA (0),
185  fD0RPResKUpgSA (0),
186  fD0RPResPiUpgSA(0),
187  fD0RPResEUpgSA(0),
188  fPt1ResPUpgSA (0),
189  fPt1ResKUpgSA (0),
190  fPt1ResPiUpgSA (0),
191  fPt1ResEUpgSA (0),
192  fRunInVertexing(isRunInVertexing),
193  fImproveTracks(kTRUE),
194  fUpdateSecVertCovMat(kTRUE),
195  fUpdateSTCovMatrix(kTRUE),
196  fUpdatePulls(kTRUE),
197  fMimicData(kFALSE),
198  fIsAOD (kTRUE),
199  fMCs (0),
200  fDebugOutput (0),
201  fDebugNtuple (0),
202  fDebugVars (0),
203  fNDebug (ndebug)
204 {
205  //
206  // Constructor to be used to create the task.
207  // The the URIs specify the resolution files to be used.
208  // They are expected to contain TGraphs with the resolutions
209  // for the current and the upgraded ITS (see code for details).
210  // One may also specify for how many tracks debug information
211  // is written to the output.
212  //
213  for(Int_t jh=0; jh<2; jh++){
214  for(Int_t ih=0; ih<4; ih++){
215  fD0RPMeanPCur[jh][ih]=0x0;
216  fD0RPMeanKCur[jh][ih]=0x0;
217  fD0RPMeanPiCur[jh][ih]=0x0;
218  fD0RPMeanECur[jh][ih]=0x0;
219  fD0RPMeanPUpg[jh][ih]=0x0;
220  fD0RPMeanKUpg[jh][ih]=0x0;
221  fD0RPMeanPiUpg[jh][ih]=0x0;
222  fD0RPMeanEUpg[jh][ih]=0x0;
223  }
224  }
225 
226  TString resfileCurURI = Form("alien:///alice/cern.ch/user/p/pwg_hf/common/Improver/%s/%s/ITSgraphs_Current.root",period,systematic);
227  TFile *resfileCur=TFile::Open(resfileCurURI.Data());
228  if(resfileCur) {
229  if(resfileCur->Get("D0RPResP" )) {
230  fD0RPResPCur =(TGraph*)(resfileCur->Get("D0RPResP" )->Clone("D0RPResPCur" ));
231  }
232  if(resfileCur->Get("D0RPResK" )) {
233  fD0RPResKCur =(TGraph*)(resfileCur->Get("D0RPResK" )->Clone("D0RPResKCur" ));
234  }
235  if(resfileCur->Get("D0RPResPi")) {
236  fD0RPResPiCur=(TGraph*)(resfileCur->Get("D0RPResPi")->Clone("D0RPResPiCur"));
237  }
238  if(resfileCur->Get("D0RPResE")) {
239  fD0RPResECur=(TGraph*)(resfileCur->Get("D0RPResE")->Clone("D0RPResECur"));
240  }
241  if(resfileCur->Get("D0RPSigmaPullRatioP" )) {
242  fD0RPSigmaPullRatioP =(TGraph*)(resfileCur->Get("D0RPSigmaPullRatioP" ));
243  }
244  if(resfileCur->Get("D0RPSigmaPullRatioK" )) {
245  fD0RPSigmaPullRatioK =(TGraph*)(resfileCur->Get("D0RPSigmaPullRatioK" ));
246  }
247  if(resfileCur->Get("D0RPSigmaPullRatioPi")) {
248  fD0RPSigmaPullRatioPi=(TGraph*)(resfileCur->Get("D0RPSigmaPullRatioPi"));
249  }
250  if(resfileCur->Get("D0RPSigmaPullRatioE")) {
251  fD0RPSigmaPullRatioE=(TGraph*)(resfileCur->Get("D0RPSigmaPullRatioE"));
252  }
253  for(Int_t j=0; j<2; j++){
254  for(Int_t i=0; i<4; i++){
255  if(resfileCur->Get(Form("D0RPMeanP_B%d_phi%d",j,i))) {
256  fD0RPMeanPCur[j][i]=(TGraph*)(resfileCur->Get(Form("D0RPMeanP_B%d_phi%d",j,i))->Clone(Form("D0RPMeanPCur_B%d_phi%d",j,i)));
257  }
258  if(resfileCur->Get(Form("D0RPMeanK_B%d_phi%d",j,i))) {
259  fD0RPMeanKCur[j][i]=(TGraph*)(resfileCur->Get(Form("D0RPMeanK_B%d_phi%d",j,i))->Clone(Form("D0RPMeanKCur_B%d_phi%d",j,i)));
260  }
261  if(resfileCur->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))) {
262  fD0RPMeanPiCur[j][i]=(TGraph*)(resfileCur->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))->Clone(Form("D0RPMeanPiCur_B%d_phi%d",j,i)));
263  }
264  if(resfileCur->Get(Form("D0RPMeanE_B%d_phi%d",j,i))) {
265  fD0RPMeanECur[j][i]=(TGraph*)(resfileCur->Get(Form("D0RPMeanE_B%d_phi%d",j,i))->Clone(Form("D0RPMeanECur_B%d_phi%d",j,i)));
266  }
267  }
268  }
269  if(resfileCur->Get("D0ZResP" )) {
270  fD0ZResPCur =(TGraph*)(resfileCur->Get("D0ZResP" )->Clone("D0ZResPCur" ));
271  }
272  if(resfileCur->Get("D0ZResK" )) {
273  fD0ZResKCur =(TGraph*)(resfileCur->Get("D0ZResK" )->Clone("D0ZResKCur" ));
274  }
275  if(resfileCur->Get("D0ZResPi" )) {
276  fD0ZResPiCur =(TGraph*)(resfileCur->Get("D0ZResPi" )->Clone("D0ZResPiCur" ));
277  }
278  if(resfileCur->Get("D0ZResE" )) {
279  fD0ZResECur =(TGraph*)(resfileCur->Get("D0ZResE" )->Clone("D0ZResECur" ));
280  }
281  if(resfileCur->Get("Pt1ResP" )) {
282  fPt1ResPCur =(TGraph*)(resfileCur->Get("Pt1ResP" )->Clone("Pt1ResPCur" ));
283  }
284  if(resfileCur->Get("Pt1ResK" )) {
285  fPt1ResKCur =(TGraph*)(resfileCur->Get("Pt1ResK" )->Clone("Pt1ResKCur" ));
286  }
287  if(resfileCur->Get("Pt1ResPi" )) {
288  fPt1ResPiCur =(TGraph*)(resfileCur->Get("Pt1ResPi" )->Clone("Pt1ResPiCur" ));
289  }
290  if(resfileCur->Get("Pt1ResE" )) {
291  fPt1ResECur =(TGraph*)(resfileCur->Get("Pt1ResE" )->Clone("Pt1ResECur" ));
292  }
293  if(resfileCur->Get("D0RPResPSA" )) {
294  fD0RPResPCurSA =(TGraph*)(resfileCur->Get("D0RPResPSA" )->Clone("D0RPResPCurSA" ));
295  }
296  if(resfileCur->Get("D0RPResKSA" )) {
297  fD0RPResKCurSA =(TGraph*)(resfileCur->Get("D0RPResKSA" )->Clone("D0RPResKCurSA" ));
298  }
299  if(resfileCur->Get("D0RPResPiSA")) {
300  fD0RPResPiCurSA=(TGraph*)(resfileCur->Get("D0RPResPiSA")->Clone("D0RPResPiCurSA"));
301  }
302  if(resfileCur->Get("D0RPResESA")) {
303  fD0RPResECurSA=(TGraph*)(resfileCur->Get("D0RPResESA")->Clone("D0RPResECurSA"));
304  }
305  if(resfileCur->Get("D0ZResPSA" )) {
306  fD0ZResPCurSA =(TGraph*)(resfileCur->Get("D0ZResPSA" )->Clone("D0ZResPCurSA" ));
307  }
308  if(resfileCur->Get("D0ZResKSA" )) {
309  fD0ZResKCurSA =(TGraph*)(resfileCur->Get("D0ZResKSA" )->Clone("D0ZResKCurSA" ));
310  }
311  if(resfileCur->Get("D0ZResPiSA" )) {
312  fD0ZResPiCurSA =(TGraph*)(resfileCur->Get("D0ZResPiSA" )->Clone("D0ZResPiCurSA" ));
313  }
314  if(resfileCur->Get("D0ZResESA" )) {
315  fD0ZResECurSA =(TGraph*)(resfileCur->Get("D0ZResESA" )->Clone("D0ZResECurSA" ));
316  }
317  if(resfileCur->Get("Pt1ResPSA" )) {
318  fPt1ResPCurSA =(TGraph*)(resfileCur->Get("Pt1ResPSA" )->Clone("Pt1ResPCurSA" ));
319  }
320  if(resfileCur->Get("Pt1ResKSA" )) {
321  fPt1ResKCurSA =(TGraph*)(resfileCur->Get("Pt1ResKSA" )->Clone("Pt1ResKCurSA" ));
322  }
323  if(resfileCur->Get("Pt1ResPiSA" )) {
324  fPt1ResPiCurSA =(TGraph*)(resfileCur->Get("Pt1ResPiSA" )->Clone("Pt1ResPiCurSA" ));
325  }
326  if(resfileCur->Get("Pt1ResESA" )) {
327  fPt1ResECurSA =(TGraph*)(resfileCur->Get("Pt1ResESA" )->Clone("Pt1ResECurSA" ));
328  }
329  delete resfileCur;
330  }
331 
332  TString resfileUpgURI = Form("alien:///alice/cern.ch/user/p/pwg_hf/common/Improver/%s/%s/ITSgraphs_NewAll-X0.3-Res4um.root",period,systematic);
333  TFile *resfileUpg=TFile::Open(resfileUpgURI.Data());
334  if(resfileUpg) {
335  if(resfileUpg->Get("D0RPResP" )) {
336  fD0RPResPUpg =(TGraph*)(resfileUpg->Get("D0RPResP" )->Clone("D0RPResPUpg" ));
337  }
338  if(resfileUpg->Get("D0RPResK" )) {
339  fD0RPResKUpg =(TGraph*)(resfileUpg->Get("D0RPResK" )->Clone("D0RPResKUpg" ));
340  }
341  if(resfileUpg->Get("D0RPResPi")) {
342  fD0RPResPiUpg=(TGraph*)(resfileUpg->Get("D0RPResPi")->Clone("D0RPResPiUpg"));
343  }
344  if(resfileUpg->Get("D0RPResE")) {
345  fD0RPResEUpg=(TGraph*)(resfileUpg->Get("D0RPResE")->Clone("D0RPResEUpg"));
346  }
347  for(Int_t j=0; j<2; j++){
348  for(Int_t i=0; i<4; i++){
349  if(resfileUpg->Get(Form("D0RPMeanP_B%d_phi%d",j,i))) {
350  fD0RPMeanPUpg[j][i]=(TGraph*)(resfileUpg->Get(Form("D0RPMeanP_B%d_phi%d",j,i))->Clone(Form("D0RPMeanPUpg_B%d_phi%d",j,i)));
351  }
352  if(resfileUpg->Get(Form("D0RPMeanK_B%d_phi%d",j,i))) {
353  fD0RPMeanKUpg[j][i]=(TGraph*)(resfileUpg->Get(Form("D0RPMeanK_B%d_phi%d",j,i))->Clone(Form("D0RPMeanKUpg_B%d_phi%d",j,i)));
354  }
355  if(resfileUpg->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))) {
356  fD0RPMeanPiUpg[j][i]=(TGraph*)(resfileUpg->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))->Clone(Form("D0RPMeanPiUpg_B%d_phi%d",j,i)));
357  }
358  if(resfileUpg->Get(Form("D0RPMeanE_B%d_phi%d",j,i))) {
359  fD0RPMeanEUpg[j][i]=(TGraph*)(resfileUpg->Get(Form("D0RPMeanE_B%d_phi%d",j,i))->Clone(Form("D0RPMeanEUpg_B%d_phi%d",j,i)));
360  }
361  }
362  }
363  if(resfileUpg->Get("D0ZResP" )) {
364  fD0ZResPUpg =(TGraph*)(resfileUpg->Get("D0ZResP" )->Clone("D0ZResPUpg" ));
365  }
366  if(resfileUpg->Get("D0ZResK" )) {
367  fD0ZResKUpg =(TGraph*)(resfileUpg->Get("D0ZResK" )->Clone("D0ZResKUpg" ));
368  }
369  if(resfileUpg->Get("D0ZResPi" )) {
370  fD0ZResPiUpg =(TGraph*)(resfileUpg->Get("D0ZResPi" )->Clone("D0ZResPiUpg" ));
371  }
372  if(resfileUpg->Get("D0ZResE" )) {
373  fD0ZResEUpg =(TGraph*)(resfileUpg->Get("D0ZResE" )->Clone("D0ZResEUpg" ));
374  }
375  if(resfileUpg->Get("Pt1ResP" )) {
376  fPt1ResPUpg =(TGraph*)(resfileUpg->Get("Pt1ResP" )->Clone("Pt1ResPUpg" ));
377  }
378  if(resfileUpg->Get("Pt1ResK" )) {
379  fPt1ResKUpg =(TGraph*)(resfileUpg->Get("Pt1ResK" )->Clone("Pt1ResKUpg" ));
380  }
381  if(resfileUpg->Get("Pt1ResPi" )) {
382  fPt1ResPiUpg =(TGraph*)(resfileUpg->Get("Pt1ResPi" )->Clone("Pt1ResPiUpg" ));
383  }
384  if(resfileUpg->Get("Pt1ResE" )) {
385  fPt1ResEUpg =(TGraph*)(resfileUpg->Get("Pt1ResE" )->Clone("Pt1ResEUpg" ));
386  }
387  if(resfileUpg->Get("D0RPResPSA" )) {
388  fD0RPResPUpgSA =(TGraph*)(resfileUpg->Get("D0RPResPSA" )->Clone("D0RPResPUpgSA" ));
389  }
390  if(resfileUpg->Get("D0RPResKSA" )) {
391  fD0RPResKUpgSA =(TGraph*)(resfileUpg->Get("D0RPResKSA" )->Clone("D0RPResKUpgSA" ));
392  }
393  if(resfileUpg->Get("D0RPResPiSA")) {
394  fD0RPResPiUpgSA=(TGraph*)(resfileUpg->Get("D0RPResPiSA")->Clone("D0RPResPiUpgSA"));
395  }
396  if(resfileUpg->Get("D0RPResESA")) {
397  fD0RPResEUpgSA=(TGraph*)(resfileUpg->Get("D0RPResESA")->Clone("D0RPResEUpgSA"));
398  }
399  if(resfileUpg->Get("D0ZResPSA" )) {
400  fD0ZResPUpgSA =(TGraph*)(resfileUpg->Get("D0ZResPSA" )->Clone("D0ZResPUpgSA" ));
401  }
402  if(resfileUpg->Get("D0ZResKSA" )) {
403  fD0ZResKUpgSA =(TGraph*)(resfileUpg->Get("D0ZResKSA" )->Clone("D0ZResKUpgSA" ));
404  }
405  if(resfileUpg->Get("D0ZResPiSA" )) {
406  fD0ZResPiUpgSA =(TGraph*)(resfileUpg->Get("D0ZResPiSA" )->Clone("D0ZResPiUpgSA" ));
407  }
408  if(resfileUpg->Get("D0ZResESA" )) {
409  fD0ZResEUpgSA =(TGraph*)(resfileUpg->Get("D0ZResESA" )->Clone("D0ZResEUpgSA" ));
410  }
411  if(resfileUpg->Get("Pt1ResPSA" )) {
412  fPt1ResPUpgSA =(TGraph*)(resfileUpg->Get("Pt1ResPSA" )->Clone("Pt1ResPUpgSA" ));
413  }
414  if(resfileUpg->Get("Pt1ResKSA" )) {
415  fPt1ResKUpgSA =(TGraph*)(resfileUpg->Get("Pt1ResKSA" )->Clone("Pt1ResKUpgSA" ));
416  }
417  if(resfileUpg->Get("Pt1ResPiSA" )) {
418  fPt1ResPiUpgSA =(TGraph*)(resfileUpg->Get("Pt1ResPiSA" )->Clone("Pt1ResPiUpgSA" ));
419  }
420  if(resfileUpg->Get("Pt1ResESA" )) {
421  fPt1ResEUpgSA =(TGraph*)(resfileUpg->Get("Pt1ResESA" )->Clone("Pt1ResEUpgSA" ));
422  }
423  delete resfileUpg;
424  }
425 
426  DefineOutput(1,TList::Class());
427 }
428 
430  //
431  // Destructor.
432  //
433  if (fDebugOutput) delete fDebugOutput;
434 }
435 
437  //
438  // Creation of user output objects.
439  //
440  fDebugOutput=new TList();
441  fDebugOutput->SetOwner();
442  fDebugOutput->SetName("debug");
443  fDebugNtuple=new TNtuple("fDebugNtuple","Smearing","pdg:ptmc:d0rpo:d0zo:pt1o:sd0rpo:sd0zo:spt1o:d0rpn:d0zn:pt1n:sd0rpn:sd0zn:spt1n:d0rpmc:d0zmc:pt1mc:pullcorr:d0zoinsigma:d0zninsigma:d0rpoinsigma:d0rpninsigma");
444  fDebugVars=new Float_t[fDebugNtuple->GetNvar()];
445 
446  fDebugOutput->Add(fDebugNtuple );
447 
456  for(Int_t j=0; j<2; j++){
457  for(Int_t i=0; i<4; i++){
458  if(fD0RPMeanPCur[j][i]) fDebugOutput->Add(fD0RPMeanPCur[j][i]);
459  if(fD0RPMeanKCur[j][i]) fDebugOutput->Add(fD0RPMeanKCur[j][i]);
460  if(fD0RPMeanPiCur[j][i]) fDebugOutput->Add(fD0RPMeanPiCur[j][i]);
461  if(fD0RPMeanECur[j][i]) fDebugOutput->Add(fD0RPMeanECur[j][i]);
462  if(fD0RPMeanPUpg[j][i]) fDebugOutput->Add(fD0RPMeanPUpg[j][i]);
463  if(fD0RPMeanKUpg[j][i]) fDebugOutput->Add(fD0RPMeanKUpg[j][i]);
464  if(fD0RPMeanPiUpg[j][i]) fDebugOutput->Add(fD0RPMeanPiUpg[j][i]);
465  if(fD0RPMeanEUpg[j][i]) fDebugOutput->Add(fD0RPMeanEUpg[j][i]);
466  }
467  }
488 
489  PostData(1,fDebugOutput);
490 }
491 
493  //
494  // The event loop
495  //
496  AliAODEvent *ev=0x0;
497  AliESDEvent *evesd=0x0;
498  Double_t bz=0.;
499 
500  if(fIsAOD) {
501  if(!fRunInVertexing) {
502  ev=dynamic_cast<AliAODEvent*>(InputEvent());
503  } else {
504  if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());
505  }
506  if(!ev) return;
507  bz=ev->GetMagneticField();
508  }
509  else {
510  evesd = dynamic_cast<AliESDEvent*>(InputEvent());
511  if (!evesd) {
512  AliError("event not found. Nothing done!");
513  return;
514  }
515  bz=evesd->GetMagneticField();
516  }
517 
518  if(fIsAOD) {
519 
520  // first loop on candidates to fill them in case of reduced AODs
521  // this is done to have the same behaviour of the improver with full (pp, p-Pb) and recuced (Pb-Pb) candidates
523 
524  // D0->Kpi
525  TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
526  if (array2Prong) {
527  for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
528  AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
529  vHF->FillRecoCand(ev,(AliAODRecoDecayHF2Prong*)decay);
530  }
531  }
532  // Dstar->Kpipi
533  TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
534  if (arrayCascade) {
535  for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
536  AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
537  vHF->FillRecoCasc(ev,((AliAODRecoCascadeHF*)decayDstar),kTRUE);
538  }
539  }
540  // Three prong
541  TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
542  if (array3Prong) {
543  for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
544  AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
545  vHF->FillRecoCand(ev,(AliAODRecoDecayHF3Prong*)decay);
546  }
547  }
548 
549 
550  // Smear all tracks
551  fMCs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
552  if (!fMCs) return;
553  if (fImproveTracks) {
554  for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
555  AliAODTrack * trk = dynamic_cast<AliAODTrack*>(ev->GetTrack(itrack));
556  if(!trk) AliFatal("Not a standard AOD");
557  SmearTrack(trk,bz);
558  }
559  }
560 
561  // TODO: recalculated primary vertex
562  AliVVertex *primaryVertex=ev->GetPrimaryVertex();
563 
564 
565  // Recalculate all candidates
566  // D0->Kpi
567  if (array2Prong) {
568  for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
569  AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
570  if(!vHF->FillRecoCand(ev,(AliAODRecoDecayHF2Prong*)decay))continue;
571 
572  // recalculate vertices
573  AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
574 
575 
576  AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
577  AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
578 
579  TObjArray ta12;
580 
581  ta12.Add(&et1); ta12.Add(&et2);
582  AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
583 
584 
585  // update secondary vertex
586  Double_t pos[3];
587  Double_t covpos[6];
588  v12->GetXYZ(pos);
589  v12->GetCovMatrix(covpos);
590  decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
591  if(fUpdateSecVertCovMat) decay->GetSecondaryVtx()->SetCovMatrix(covpos);
592  decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
593 
594  // update d0
595  Double_t d0z0[2],covd0z0[3];
596  Double_t d0[2],d0err[2];
597  et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
598  d0[0]=d0z0[0];
599  d0err[0] = TMath::Sqrt(covd0z0[0]);
600  et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
601  d0[1]=d0z0[0];
602  d0err[1] = TMath::Sqrt(covd0z0[0]);
603  decay->Setd0Prongs(2,d0);
604  decay->Setd0errProngs(2,d0err);
605  //
606 
607 
608  Double_t xdummy=0.,ydummy=0.;
609  Double_t dca;
610  dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
611  decay->SetDCA(dca);
612 
613 
614 
615  Double_t px[2],py[2],pz[2];
616  for (Int_t i=0;i<2;++i) {
617  AliExternalTrackParam et;
618  et.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(i)));
619  et.PropagateToDCA(v12,bz,100.,d0z0,covd0z0);
620  px[i]=et.Px();
621  py[i]=et.Py();
622  pz[i]=et.Pz();
623  }
624  decay->SetPxPyPzProngs(2,px,py,pz);
625  delete v12;
626  }
627  }
628 
629 
630  // Dstar->Kpipi
631  if (arrayCascade) {
632  for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
633  AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
634  if(!vHF->FillRecoCasc(ev,((AliAODRecoCascadeHF*)decayDstar),kTRUE))continue;
635  //Get D0 from D*
637 
638  // recalculate vertices
639  //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
640 
641  //soft pion
642  AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
643 
644  //track D0
645  AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
646 
648 
649  // update d0
650  Double_t d0z0[2],covd0z0[3];
651  Double_t d01[2],d01err[2];
652 
653  //the D*
654  et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
655  d01[0]=d0z0[0];
656  d01err[0] = TMath::Sqrt(covd0z0[0]);
657  trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
658  d01[1]=d0z0[0];
659  d01err[1] = TMath::Sqrt(covd0z0[0]);
660  decayDstar->Setd0Prongs(2,d01);
661  decayDstar->Setd0errProngs(2,d01err);
662 
663  // delete v12;
664  delete trackD0; trackD0=NULL;
665 
666  // a run for D*
667  Double_t px1[2],py1[2],pz1[2];
668  for (Int_t i=0;i<2;++i) {
669  const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));
670  px1[i]=t1->Px();
671  py1[i]=t1->Py();
672  pz1[i]=t1->Pz();
673  }
674  decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
675 
676  }
677  }
678 
679 
680  // Three prong
681  if (array3Prong) {
682  for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
683  AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
684  if(!vHF->FillRecoCand(ev,(AliAODRecoDecayHF3Prong*)decay))continue;
685 
686  // recalculate vertices
687  AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
688  AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
689  AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
690  AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
691  TObjArray ta123,ta12,ta23;
692  ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
693  ta12. Add(&et1);ta12 .Add(&et2);
694  ta23 .Add(&et2);ta23 .Add(&et3);
695  AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
696  AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
697  AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
698 
699  // update secondary vertex
700  Double_t pos[3];
701  Double_t covpos[6];
702  v123->GetXYZ(pos);
703  v123->GetCovMatrix(covpos);
704  decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
705  if(fUpdateSecVertCovMat) decay->GetSecondaryVtx()->SetCovMatrix(covpos);
706  decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
707 
708  // update d0 for all progs
709  Double_t d0z0[2],covd0z0[3];
710  Double_t d0[3],d0err[3];
711  et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
712  d0[0]=d0z0[0];
713  d0err[0] = TMath::Sqrt(covd0z0[0]);
714  et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
715  d0[1]=d0z0[0];
716  d0err[1] = TMath::Sqrt(covd0z0[0]);
717  et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
718  d0[2]=d0z0[0];
719  d0err[2] = TMath::Sqrt(covd0z0[0]);
720  decay->Setd0Prongs (3,d0 );
721  decay->Setd0errProngs(3,d0err);
722  // TODO: setter missing
723 
724  // update dca for prong combinations
725  Double_t xdummy=0.,ydummy=0.;
726  Double_t dca[3];
727  dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
728  dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
729  dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
730  decay->SetDCAs(3,dca);
731 
732  // update sigmavertex = dispersion
733  Float_t sigmaV=v123->GetDispersion();
734  decay->SetSigmaVert(sigmaV);
735  // update dist12 and dist23
736  primaryVertex->GetXYZ(pos);
737  decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
738  +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
739  +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
740  decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
741  +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
742  +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
743 
744 
745  Double_t px[3],py[3],pz[3];
746  for (Int_t i=0;i<3;++i) {
747  AliExternalTrackParam et;
748  et.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(i)));
749  et.PropagateToDCA(v123,bz,100.,d0z0,covd0z0);
750  px[i]=et.Px();
751  py[i]=et.Py();
752  pz[i]=et.Pz();
753  }
754  decay->SetPxPyPzProngs(3,px,py,pz);
755 
756  delete v123;delete v12;delete v23;
757  }
758  }
759  delete vHF;
760 
761  } // end AOD
762  else {
763  //
764  // In case of ESD: only smear all tracks
765  //
766  if (!fMCEvent) return;
767  if (fImproveTracks) {
768  for(Int_t itrack=0;itrack<evesd->GetNumberOfTracks();++itrack) {
769  AliESDtrack * trk = dynamic_cast<AliESDtrack*>(evesd->GetTrack(itrack));
770  if(!trk) AliFatal("No a standard ESD");
771  SmearTrack(trk,bz);
772  }
773  }
774  }// end ESD
775 
776 }
777 
779  // Early exit, if this track has nothing in common with the ITS
780 
781  if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
782  return;
783 
784  // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))
785  if (TESTBIT(track->GetITSClusterMap(),7)) return;
786  //
787 
788 
789  // Get reconstructed track parameters
790  AliExternalTrackParam et; et.CopyFromVTrack(track);
791  Double_t *param=const_cast<Double_t*>(et.GetParameter());
792 
793 
794  Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
795 
796  // Get MC info
797  Int_t imc=track->GetLabel();
798  if (imc<=0) return;
799  Double_t mcx[3];
800  Double_t mcp[3];
801  Double_t mccv[36]={0.};
802  Short_t mcc;
803  const AliVParticle *mc= 0x0;
804  if(fIsAOD) {
805  mc = dynamic_cast<AliVParticle*>(fMCs->At(imc));
806  }
807  else {
808  mc = dynamic_cast<AliVParticle*>(fMCEvent->GetTrack(imc));
809  }
810  if(!mc) return;
811  mc->XvYvZv(mcx);
812  mc->PxPyPz(mcp);
813  mcc=mc->Charge();
814  AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
815  const Double_t *parammc=mct.GetParameter();
816 //TODO: const Double_t *covermc=mct.GetCovariance();
817  AliVertex vtx(mcx,1.,1);
818 
819  // Correct reference points and frames according to MC
820  // TODO: B-Field correct?
821  // TODO: failing propagation....
822  et.PropagateToDCA(&vtx,track->GetBz(),10.);
823  et.Rotate(mct.GetAlpha());
824 
825  // Select appropriate smearing functions
826  Double_t ptmc=TMath::Abs(mc->Pt());
827  Double_t phimc=mc->Phi();
828  Int_t phiBin=PhiBin(phimc);
829  Int_t magfield=0;
830  if(bz<0.) magfield=0;
831  else if(bz>0.)magfield=1;
832  Double_t sd0rpn=0.;
833  Double_t sd0mrpn=0.;
834  Double_t sd0zn =0.;
835  Double_t spt1n =0.;
836  Double_t sd0rpo=0.;
837  Double_t sd0mrpo=0.;
838  Double_t sd0zo =0.;
839  Double_t spt1o =0.;
840  Double_t pullcorr=1.;
841  Int_t pdgcode = 0;
842  if(fIsAOD){
843  const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mc);
844  if(!amcpart) return;
845  pdgcode = amcpart->GetPdgCode();
846  } else {
847  const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mc);
848  if(!emcpart) return;
849  pdgcode = emcpart->PdgCode();
850  }
851 
852  switch (pdgcode) {
853  case 2212: case -2212:
855  sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);
856  spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);
858  sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);
859  spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);
860  sd0mrpo=EvalGraph(ptmc,fD0RPMeanPCur[magfield][phiBin],fD0RPMeanPCur[magfield][phiBin]);
861  sd0mrpn=EvalGraph(ptmc,fD0RPMeanPUpg[magfield][phiBin],fD0RPMeanPUpg[magfield][phiBin]);
863  break;
864  case 321: case -321:
866  sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);
867  spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);
869  sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);
870  spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);
871  sd0mrpo=EvalGraph(ptmc,fD0RPMeanKCur[magfield][phiBin],fD0RPMeanKCur[magfield][phiBin]);
872  sd0mrpn=EvalGraph(ptmc,fD0RPMeanKUpg[magfield][phiBin],fD0RPMeanKUpg[magfield][phiBin]);
874  break;
875  case 211: case -211:
882  sd0mrpo=EvalGraph(ptmc,fD0RPMeanPiCur[magfield][phiBin],fD0RPMeanPiCur[magfield][phiBin]);
883  sd0mrpn=EvalGraph(ptmc,fD0RPMeanPiUpg[magfield][phiBin],fD0RPMeanPiUpg[magfield][phiBin]);
885  break;
886  case 11: case -11:
888  sd0zo =EvalGraph(ptmc,fD0ZResECur,fD0ZResECurSA);
889  spt1o =EvalGraph(ptmc,fPt1ResECur,fPt1ResECurSA);
891  sd0zn =EvalGraph(ptmc,fD0ZResEUpg,fD0ZResEUpgSA);
892  spt1n =EvalGraph(ptmc,fPt1ResEUpg,fPt1ResEUpgSA);
893  sd0mrpo=EvalGraph(ptmc,fD0RPMeanECur[magfield][phiBin],fD0RPMeanECur[magfield][phiBin]);
894  sd0mrpn=EvalGraph(ptmc,fD0RPMeanEUpg[magfield][phiBin],fD0RPMeanEUpg[magfield][phiBin]);
896  break;
897  default:
898  return;
899  }
900 
901  // Use the same units (i.e. cm and GeV/c)! TODO: pt!
902  sd0rpo*=1.e-4;
903  sd0zo *=1.e-4;
904  sd0rpn*=1.e-4;
905  sd0zn *=1.e-4;
906  sd0mrpo*=1.e-4;
907  sd0mrpn*=1.e-4;
908 
909  // Apply the smearing
910  Double_t d0zo =param [1];
911  Double_t d0zmc =parammc[1];
912  Double_t d0rpo =param [0];
913  Double_t d0rpmc=parammc[0];
914  Double_t pt1o =param [4];
915  Double_t pt1mc =parammc[4];
916  Double_t dd0zo =d0zo-d0zmc;
917  Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
918  Double_t d0zn =d0zmc+dd0zn;
919  Double_t dd0rpo=d0rpo-d0rpmc;
920  Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
921  Double_t dd0mrpn=TMath::Abs(sd0mrpn)-TMath::Abs(sd0mrpo);
922  Double_t d0rpn =d0rpmc+dd0rpn-dd0mrpn;
923  Double_t d0zoinsigma = 0.;
924  if(covar[0] > 0.) d0zoinsigma = d0zo/TMath::Sqrt(covar[2]);
925  Double_t d0rpoinsigma = 0.;
926  if(covar[2] > 0.) d0rpoinsigma = d0rpo/TMath::Sqrt(covar[0]);
927 
928  if(fMimicData){
929  dd0mrpn=sd0mrpn-sd0mrpo;
930  d0rpn =d0rpmc+dd0rpn+dd0mrpn;
931  }
932 
933  Double_t dpt1o =pt1o-pt1mc;
934  Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
935  Double_t pt1n =pt1mc+dpt1n;
936  param[0]=d0rpn;
937  param[1]=d0zn ;
938  param[4]=pt1n ;
939 
940  //cov matrix update
941  if(fUpdateSTCovMatrix){
942  if(sd0rpo>0.) covar[0]*=(sd0rpn/sd0rpo)*(sd0rpn/sd0rpo);//yy
943  if(sd0zo>0. && sd0rpo>0.)covar[1]*=(sd0rpn/sd0rpo)*(sd0zn/sd0zo);//yz
944  if(sd0zo>0.) covar[2]*=(sd0zn/sd0zo)*(sd0zn/sd0zo);//zz
945  if(sd0rpo>0.) covar[3]*=(sd0rpn/sd0rpo);//yl
946  if(sd0zo>0.) covar[4]*=(sd0zn/sd0zo);//zl
947  if(sd0rpo>0.) covar[6]*=(sd0rpn/sd0rpo);//ysenT
948  if(sd0zo>0.) covar[7]*=(sd0zn/sd0zo);//zsenT
949  if(sd0rpo>0. && spt1o>0.)covar[10]*=(sd0rpn/sd0rpo)*(spt1n/spt1o);//ypt
950  if(sd0zo>0. && spt1o>0.) covar[11]*=(sd0zn/sd0zo)*(spt1n/spt1o);//zpt
951  if(spt1o>0.) covar[12]*=(spt1n/spt1o);//sinPhipt
952  if(spt1o>0.) covar[13]*=(spt1n/spt1o);//tanTpt
953  if(spt1o>0.) covar[14]*=(spt1n/spt1o)*(spt1n/spt1o);//ptpt
954  }
955  if(fUpdatePulls){
956 
957  covar[0]*=pullcorr*pullcorr;//yy
958  covar[1]*=pullcorr;//yz
959  covar[3]*=pullcorr;//yl
960  covar[6]*=pullcorr;//ysenT
961  covar[10]*=pullcorr;//ypt
962 
963  }
964 
965  Double_t d0zninsigma = 0.;
966  if(covar[0] > 0.) d0zninsigma = d0zn/TMath::Sqrt(covar[2]);
967  Double_t d0rpninsigma = 0.;
968  if(covar[2] > 0.) d0rpninsigma = d0rpn/TMath::Sqrt(covar[0]);
969 
970  // Copy the smeared parameters to the AOD track
971  Double_t x[3];
972  Double_t p[3];
973  et.GetXYZ(x);
974  et.GetPxPyPz(p);
975  Double_t cv[21];
976  et.GetCovarianceXYZPxPyPz(cv);
977  if(fIsAOD) {
978  AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
979  aodtrack->SetPosition(x,kFALSE);
980  aodtrack->SetP(p,kTRUE);
981  aodtrack->SetCovMatrix(cv);
982  // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
983  UChar_t itsClusterMap = aodtrack->GetITSClusterMap();
984  SETBIT(itsClusterMap,7);
985  aodtrack->SetITSClusterMap(itsClusterMap);
986  //
987 
988  } else {
989  AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
990  Short_t sign = esdtrack->Charge();
991  esdtrack->Set(x,p,cv,sign);
992  esdtrack->RelateToVVertex(InputEvent()->GetPrimaryVertex(), bz,100.);
993  // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
994  UChar_t itsClusterMap = esdtrack->GetITSClusterMap();
995  SETBIT(itsClusterMap,7);
996  esdtrack->SetITSClusterMap(itsClusterMap);
997  }
998 
999 
1000  // write out debug infos
1001  if (fDebugNtuple->GetEntries()<fNDebug) {
1002  Int_t idbg=0;
1003  fDebugVars[idbg++]=pdgcode;
1004  fDebugVars[idbg++]=ptmc ;
1005  fDebugVars[idbg++]=d0rpo ;
1006  fDebugVars[idbg++]=d0zo ;
1007  fDebugVars[idbg++]=pt1o ;
1008  fDebugVars[idbg++]=sd0rpo;
1009  fDebugVars[idbg++]=sd0zo ;
1010  fDebugVars[idbg++]=spt1o ;
1011  fDebugVars[idbg++]=d0rpn ;
1012  fDebugVars[idbg++]=d0zn ;
1013  fDebugVars[idbg++]=pt1n ;
1014  fDebugVars[idbg++]=sd0rpn;
1015  fDebugVars[idbg++]=sd0zn ;
1016  fDebugVars[idbg++]=spt1n ;
1017  fDebugVars[idbg++]=d0rpmc;
1018  fDebugVars[idbg++]=d0zmc ;
1019  fDebugVars[idbg++]=pt1mc ;
1020  fDebugVars[idbg++]=pullcorr ;
1021  fDebugVars[idbg++]=d0zoinsigma ;
1022  fDebugVars[idbg++]=d0zninsigma ;
1023  fDebugVars[idbg++]=d0rpoinsigma ;
1024  fDebugVars[idbg++]=d0rpninsigma ;
1025  fDebugNtuple->Fill(fDebugVars);
1026  PostData(1,fDebugOutput);
1027  }
1028 }
1029 
1030 AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
1031  //
1032  // Helper function to recalculate a vertex.
1033  //
1034 
1035  static UShort_t ids[]={1,2,3}; //TODO: unsave...
1036  AliVertexerTracks vertexer(bField);
1037  vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
1038  AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
1039  return vertex;
1040 }
1041 
1043  //
1044  // Evaluates a TGraph without linear extrapolation. Instead the last
1045  // valid point of the graph is used when out of range.
1046  // The function assumes an ascending order of X.
1047  //
1048 
1049  if(!graph) return 0.;
1050 
1051  // TODO: find a pretty solution for this:
1052  Int_t n =graph->GetN();
1053  Double_t xmin=graph->GetX()[0 ];
1054  Double_t xmax=graph->GetX()[n-1];
1055  if (x<xmin) {
1056  if(!graphSA) return graph->Eval(xmin);
1057  Double_t xminSA=graphSA->GetX()[0];
1058  if(x<xminSA) return graphSA->Eval(xminSA);
1059  return graphSA->Eval(x);
1060  }
1061  if (x>xmax) return graph->Eval(xmax);
1062  return graph->Eval(x);
1063 }
1064 
1065 //________________________________________________________________________
1066 
1068  Double_t pi=TMath::Pi();
1069  if(phi>2.*pi || phi<0.) return -1;
1070  if((phi<=(pi/4.)) || (phi>7.*(pi/4.))) return 0;
1071  if((phi>(pi/4.)) && (phi<=3.*(pi/4.))) return 1;
1072  if((phi>3.*(pi/4.)) && (phi<=5.*(pi/4.))) return 2;
1073  if((phi>(5.*pi/4.)) && (phi<=7.*(pi/4.))) return 3;
1074  return -1;
1075 }
1076 
1077 ClassImp(AliAnalysisTaskSEImproveITS);
1078 
TGraph * fD0RPMeanEUpg[2][4]
new pt dep. d0 mean in rphi for pions in 4 phi regions
TGraph * fPt1ResPiCur
old pt dep. 1/pt res. for kaons
TGraph * fD0RPMeanPUpg[2][4]
new pt dep. d0 res. in rphi for electrons
TGraph * fD0ZResPUpg
old pt dep. 1/pt res. for electrons
TGraph * fD0RPResKCur
old pt dep. d0 res. in rphi for protons
Bool_t fMimicData
flag to switch on/off the correction of the pulls
TGraph * fD0RPResPiUpgSA
new standalone pt dep. d0 res. in rphi for kaons
double Double_t
Definition: External.C:58
TGraph * fD0ZResKUpgSA
new standalone pt dep. d0 res. in z for protons
TGraph * fD0RPResPCurSA
old standalone pt dep. d0 res. in z for electrons
virtual void UserCreateOutputObjects()
Implementation of interface methods.
TGraph * fD0RPResPiCur
old pt dep. d0 res. in rphi for kaons
TGraph * fPt1ResKUpgSA
new standalone pt dep. 1/pt res. for protons
TGraph * fD0ZResPUpgSA
old standalone pt dep. 1/pt res. for electrons
TGraph * fD0RPMeanPiCur[2][4]
old pt dep. d0 mean. in rphi for kaons in 4 phi regions
Double_t EvalGraph(Double_t x, const TGraph *graph, const TGraph *graphSA=0) const
Helper functions.
TGraph * fD0RPMeanPiUpg[2][4]
new pt dep. d0 mean in rphi for kaons in 4 phi regions
Double_t bz
Bool_t fUpdateSecVertCovMat
this is always kTRUE. kFALSE only if re-running on already improved AODs
TGraph * fD0ZResPiCurSA
old standalone pt dep. d0 res. in z for kaons
virtual void UserExec(Option_t *option)
TGraph * fD0RPMeanPCur[2][4]
pt dep. d0 sigma pull MC/data in rphi for electrons
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
TList * fDebugOutput
! collection of debug output
TGraph * fPt1ResPiUpg
new pt dep. 1/pt res. for kaons
TGraph * fPt1ResPCurSA
old standalone pt dep. d0 res. in rphi for electrons
TGraph * fPt1ResPiCurSA
old standalone pt dep. 1/pt res. for kaons
TGraph * fD0RPSigmaPullRatioE
pt dep. d0 sigma pull MC/data in rphi for pions
TGraph * fD0RPResECur
old pt dep. d0 res. in rphi for pions
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
TGraph * fPt1ResPUpg
new pt dep. d0 mean in rphi for electrons in 4 phi regions
TGraph * fD0RPResKCurSA
old standalone pt dep. d0 res. in rphi for protons
TGraph * fPt1ResECurSA
old standalone pt dep. 1/pt res. for pions
TGraph * fPt1ResKUpg
new pt dep. 1/pt res. for protons
int Int_t
Definition: External.C:63
TGraph * fPt1ResKCur
old pt dep. 1/pt res. for protons
TGraph * fD0RPResKUpgSA
new standalone pt dep. d0 res. in rphi for protons
float Float_t
Definition: External.C:68
TGraph * fPt1ResPUpgSA
new standalone pt dep. d0 res. in rphi for electrons
TGraph * fPt1ResECur
old pt dep. 1/pt res. for pions
AliAODTrack * GetBachelor() const
TGraph * fD0RPSigmaPullRatioK
pt dep. d0 sigma pull MC/data in rphi for protons
TGraph * fD0RPResPiCurSA
old standalone pt dep. d0 res. in rphi for kaons
TGraph * fD0RPSigmaPullRatioP
old pt dep. d0 res. in rphi for electrons
TGraph * fD0ZResKCurSA
old standalone pt dep. d0 res. in z for protons
TGraph * fD0ZResPCurSA
new pt dep. 1/pt res. for electrons
void SmearTrack(AliVTrack *track, Double_t bz)
Float_t * fDebugVars
! variables to store as degug info
TGraph * fD0RPResPCur
old pt dep. d0 res. in z for electrons
short Short_t
Definition: External.C:23
Bool_t fRunInVertexing
new standalone pt dep. 1/pt res. for electrons
TGraph * fD0ZResEUpg
new pt dep. d0 res. in z for pions
TGraph * fPt1ResPiUpgSA
new standalone pt dep. 1/pt res. for kaons
TGraph * fPt1ResEUpgSA
new standalone pt dep. 1/pt res. for pions
TGraph * fD0RPResPUpgSA
new standalone pt dep. d0 res. in z for electrons
TGraph * fPt1ResKCurSA
old standalone pt dep. 1/pt res. for protons
decay
Definition: HFPtSpectrum.C:41
TGraph * fD0ZResECurSA
old standalone pt dep. d0 res. in z for pions
TGraph * fPt1ResEUpg
new pt dep. 1/pt res. for pions
TGraph * fD0RPSigmaPullRatioPi
pt dep. d0 sigma pull MC/data in rphi for kaons
Bool_t fImproveTracks
flag to run hybrid task before the vertexingHF task or in standard mode
TGraph * fD0RPResECurSA
old standalone pt dep. d0 res. in rphi for pions
TGraph * fD0ZResKCur
old pt dep. d0 res. in z for protons
TGraph * fD0ZResEUpgSA
new standalone pt dep. d0 res. in z for pions
void Setd0errProngs(Int_t nprongs, Double_t *d0)
TGraph * fPt1ResPCur
old pt dep. d0 mean. in rphi for electrons in 4 phi regions
TGraph * fD0RPResPUpg
new pt dep. d0 res. in z for electrons
TGraph * fD0ZResECur
old pt dep. d0 res. in z for pions
Bool_t fUpdateSTCovMatrix
flag to swicth on/off the modification of the sec vert cov matrix
TGraph * fD0ZResPiUpg
new pt dep. d0 res. in z for kaons
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
TGraph * fD0ZResPiCur
old pt dep. d0 res. in z for kaons
Bool_t fUpdatePulls
flag to switch on/off the update of the single track covariance matrix
const Double_t pi
TGraph * fD0RPResEUpgSA
new standalone pt dep. d0 res. in rphi for pions
bool Bool_t
Definition: External.C:53
void SetSigmaVert(Double_t sigmaVert)
TGraph * fD0RPResEUpg
new pt dep. d0 res. in rphi for pions
AliAODRecoDecayHF2Prong * Get2Prong() const
TGraph * fD0ZResPiUpgSA
new standalone pt dep. d0 res. in z for kaons
TGraph * fD0RPResPiUpg
new pt dep. d0 res. in rphi for kaons
TGraph * fD0RPMeanKUpg[2][4]
new pt dep. d0 mean in rphi for protons in 4 phi regions
AliESDVertex * RecalculateVertex(const AliVVertex *old, TObjArray *tracks, Double_t bField)
TGraph * fD0RPMeanKCur[2][4]
old pt dep. d0 mean. in rphi for protons in 4 phi regions
TNtuple * fDebugNtuple
! debug send on output slot 1
TGraph * fD0RPMeanECur[2][4]
old pt dep. d0 mean. in rphi for pions in 4 phi regions
TGraph * fD0RPResKUpg
new pt dep. d0 res. in rphi for protons
TGraph * fD0ZResKUpg
new pt dep. d0 res. in z for protons