AliPhysics  b76e98e (b76e98e)
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 =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
231  fD0RPResPCur ->SetName("D0RPResPCur" );
232  }
233  if(resfileCur->Get("D0RPResK" )) {
234  fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
235  fD0RPResKCur ->SetName("D0RPResKCur" );
236  }
237  if(resfileCur->Get("D0RPResPi")) {
238  fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
239  fD0RPResPiCur->SetName("D0RPResPiCur");
240  }
241  if(resfileCur->Get("D0RPResE")) {
242  fD0RPResECur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResE")));
243  fD0RPResECur->SetName("D0RPResECur");
244  }
245  if(resfileCur->Get("D0RPSigmaPullRatioP" )) {
246  fD0RPSigmaPullRatioP =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPSigmaPullRatioP" )));
247  }
248  if(resfileCur->Get("D0RPSigmaPullRatioK" )) {
249  fD0RPSigmaPullRatioK =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPSigmaPullRatioK" )));
250  }
251  if(resfileCur->Get("D0RPSigmaPullRatioPi")) {
252  fD0RPSigmaPullRatioPi=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPSigmaPullRatioPi")));
253  }
254  if(resfileCur->Get("D0RPSigmaPullRatioE")) {
255  fD0RPSigmaPullRatioE=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPSigmaPullRatioE")));
256  }
257  for(Int_t j=0; j<2; j++){
258  for(Int_t i=0; i<4; i++){
259  if(resfileCur->Get(Form("D0RPMeanP_B%d_phi%d",j,i))) {
260  fD0RPMeanPCur[j][i]=new TGraph(*static_cast<TGraph*>(resfileCur->Get(Form("D0RPMeanP_B%d_phi%d",j,i))));
261  fD0RPMeanPCur[j][i]->SetName(Form("D0RPMeanPCur_B%d_phi%d",j,i));
262  }
263  if(resfileCur->Get(Form("D0RPMeanK_B%d_phi%d",j,i))) {
264  fD0RPMeanKCur[j][i]=new TGraph(*static_cast<TGraph*>(resfileCur->Get(Form("D0RPMeanK_B%d_phi%d",j,i))));
265  fD0RPMeanKCur[j][i]->SetName(Form("D0RPMeanKCur_B%d_phi%d",j,i));
266  }
267  if(resfileCur->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))) {
268  fD0RPMeanPiCur[j][i]=new TGraph(*static_cast<TGraph*>(resfileCur->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))));
269  fD0RPMeanPiCur[j][i]->SetName(Form("D0RPMeanPiCur_B%d_phi%d",j,i));
270  }
271  if(resfileCur->Get(Form("D0RPMeanE_B%d_phi%d",j,i))) {
272  fD0RPMeanECur[j][i]=new TGraph(*static_cast<TGraph*>(resfileCur->Get(Form("D0RPMeanE_B%d_phi%d",j,i))));
273  fD0RPMeanECur[j][i]->SetName(Form("D0RPMeanECur_B%d_phi%d",j,i));
274  }
275  }
276  }
277  if(resfileCur->Get("D0ZResP" )) {
278  fD0ZResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP" )));
279  fD0ZResPCur ->SetName("D0ZResPCur" );
280  }
281  if(resfileCur->Get("D0ZResK" )) {
282  fD0ZResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK" )));
283  fD0ZResKCur ->SetName("D0ZResKCur" );
284  }
285  if(resfileCur->Get("D0ZResPi" )) {
286  fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
287  fD0ZResPiCur ->SetName("D0ZResPiCur" );
288  }
289  if(resfileCur->Get("D0ZResE" )) {
290  fD0ZResECur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResE" )));
291  fD0ZResECur ->SetName("D0ZResECur" );
292  }
293  if(resfileCur->Get("Pt1ResP" )) {
294  fPt1ResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP" )));
295  fPt1ResPCur ->SetName("Pt1ResPCur" );
296  }
297  if(resfileCur->Get("Pt1ResK" )) {
298  fPt1ResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK" )));
299  fPt1ResKCur ->SetName("Pt1ResKCur" );
300  }
301  if(resfileCur->Get("Pt1ResPi" )) {
302  fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
303  fPt1ResPiCur ->SetName("Pt1ResPiCur" );
304  }
305  if(resfileCur->Get("Pt1ResE" )) {
306  fPt1ResECur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResE" )));
307  fPt1ResECur ->SetName("Pt1ResECur" );
308  }
309  if(resfileCur->Get("D0RPResPSA" )) {
310  fD0RPResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPSA" )));
311  fD0RPResPCurSA ->SetName("D0RPResPCurSA" );
312  }
313  if(resfileCur->Get("D0RPResKSA" )) {
314  fD0RPResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResKSA" )));
315  fD0RPResKCurSA ->SetName("D0RPResKCurSA" );
316  }
317  if(resfileCur->Get("D0RPResPiSA")) {
318  fD0RPResPiCurSA=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPiSA")));
319  fD0RPResPiCurSA->SetName("D0RPResPiCurSA");
320  }
321  if(resfileCur->Get("D0RPResESA")) {
322  fD0RPResECurSA=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResESA")));
323  fD0RPResECurSA->SetName("D0RPResECurSA");
324  }
325  if(resfileCur->Get("D0ZResPSA" )) {
326  fD0ZResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPSA" )));
327  fD0ZResPCurSA ->SetName("D0ZResPCurSA" );
328  }
329  if(resfileCur->Get("D0ZResKSA" )) {
330  fD0ZResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResKSA" )));
331  fD0ZResKCurSA ->SetName("D0ZResKCurSA" );
332  }
333  if(resfileCur->Get("D0ZResPiSA" )) {
334  fD0ZResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPiSA" )));
335  fD0ZResPiCurSA ->SetName("D0ZResPiCurSA" );
336  }
337  if(resfileCur->Get("D0ZResESA" )) {
338  fD0ZResECurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResESA" )));
339  fD0ZResECurSA ->SetName("D0ZResECurSA" );
340  }
341  if(resfileCur->Get("Pt1ResPSA" )) {
342  fPt1ResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPSA" )));
343  fPt1ResPCurSA ->SetName("Pt1ResPCurSA" );
344  }
345  if(resfileCur->Get("Pt1ResKSA" )) {
346  fPt1ResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResKSA" )));
347  fPt1ResKCurSA ->SetName("Pt1ResKCurSA" );
348  }
349  if(resfileCur->Get("Pt1ResPiSA" )) {
350  fPt1ResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPiSA" )));
351  fPt1ResPiCurSA ->SetName("Pt1ResPiCurSA" );
352  }
353  if(resfileCur->Get("Pt1ResESA" )) {
354  fPt1ResECurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResESA" )));
355  fPt1ResECurSA ->SetName("Pt1ResECurSA" );
356  }
357  delete resfileCur;
358  }
359 
360  TString resfileUpgURI = Form("alien:///alice/cern.ch/user/p/pwg_hf/common/Improver/%s/%s/ITSgraphs_NewAll-X0.3-Res4um.root",period,systematic);
361  TFile *resfileUpg=TFile::Open(resfileUpgURI.Data());
362  if(resfileUpg) {
363  if(resfileUpg->Get("D0RPResP" )) {
364  fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
365  fD0RPResPUpg ->SetName("D0RPResPUpg" );
366  }
367  if(resfileUpg->Get("D0RPResK" )) {
368  fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
369  fD0RPResKUpg ->SetName("D0RPResKUpg" );
370  }
371  if(resfileUpg->Get("D0RPResPi")) {
372  fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
373  fD0RPResPiUpg->SetName("D0RPResPiUpg");
374  }
375  if(resfileUpg->Get("D0RPResE")) {
376  fD0RPResEUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResE")));
377  fD0RPResEUpg->SetName("D0RPResEUpg");
378  }
379  for(Int_t j=0; j<2; j++){
380  for(Int_t i=0; i<4; i++){
381  if(resfileUpg->Get(Form("D0RPMeanP_B%d_phi%d",j,i))) {
382  fD0RPMeanPUpg[j][i]=new TGraph(*static_cast<TGraph*>(resfileUpg->Get(Form("D0RPMeanP_B%d_phi%d",j,i))));
383  fD0RPMeanPUpg[j][i]->SetName(Form("D0RPMeanPUpg_B%d_phi%d",j,i));
384  }
385  if(resfileUpg->Get(Form("D0RPMeanK_B%d_phi%d",j,i))) {
386  fD0RPMeanKUpg[j][i]=new TGraph(*static_cast<TGraph*>(resfileUpg->Get(Form("D0RPMeanK_B%d_phi%d",j,i))));
387  fD0RPMeanKUpg[j][i]->SetName(Form("D0RPMeanKUpg_B%d_phi%d",j,i));
388  }
389  if(resfileUpg->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))) {
390  fD0RPMeanPiUpg[j][i]=new TGraph(*static_cast<TGraph*>(resfileUpg->Get(Form("D0RPMeanPi_B%d_phi%d",j,i))));
391  fD0RPMeanPiUpg[j][i]->SetName(Form("D0RPMeanPiUpg_B%d_phi%d",j,i));
392  }
393  if(resfileUpg->Get(Form("D0RPMeanE_B%d_phi%d",j,i))) {
394  fD0RPMeanEUpg[j][i]=new TGraph(*static_cast<TGraph*>(resfileUpg->Get(Form("D0RPMeanE_B%d_phi%d",j,i))));
395  fD0RPMeanEUpg[j][i]->SetName(Form("D0RPMeanEUpg_B%d_phi%d",j,i));
396  }
397  }
398  }
399  if(resfileUpg->Get("D0ZResP" )) {
400  fD0ZResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP" )));
401  fD0ZResPUpg ->SetName("D0ZResPUpg" );
402  }
403  if(resfileUpg->Get("D0ZResK" )) {
404  fD0ZResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK" )));
405  fD0ZResKUpg ->SetName("D0ZResKUpg" );
406  }
407  if(resfileUpg->Get("D0ZResPi" )) {
408  fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
409  fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
410  }
411  if(resfileUpg->Get("D0ZResE" )) {
412  fD0ZResEUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResE" )));
413  fD0ZResEUpg ->SetName("D0ZResEUpg" );
414  }
415  if(resfileUpg->Get("Pt1ResP" )) {
416  fPt1ResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP" )));
417  fPt1ResPUpg ->SetName("Pt1ResPUpg" );
418  }
419  if(resfileUpg->Get("Pt1ResK" )) {
420  fPt1ResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK" )));
421  fPt1ResKUpg ->SetName("Pt1ResKUpg" );
422  }
423  if(resfileUpg->Get("Pt1ResPi" )) {
424  fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
425  fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
426  }
427  if(resfileUpg->Get("Pt1ResE" )) {
428  fPt1ResEUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResE" )));
429  fPt1ResEUpg ->SetName("Pt1ResEUpg" );
430  }
431  if(resfileUpg->Get("D0RPResPSA" )) {
432  fD0RPResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPSA" )));
433  fD0RPResPUpgSA ->SetName("D0RPResPUpgSA" );
434  }
435  if(resfileUpg->Get("D0RPResKSA" )) {
436  fD0RPResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResKSA" )));
437  fD0RPResKUpgSA ->SetName("D0RPResKUpgSA" );
438  }
439  if(resfileUpg->Get("D0RPResPiSA")) {
440  fD0RPResPiUpgSA=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPiSA")));
441  fD0RPResPiUpgSA->SetName("D0RPResPiUpgSA");
442  }
443  if(resfileUpg->Get("D0RPResESA")) {
444  fD0RPResEUpgSA=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResESA")));
445  fD0RPResEUpgSA->SetName("D0RPResEUpgSA");
446  }
447  if(resfileUpg->Get("D0ZResPSA" )) {
448  fD0ZResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPSA" )));
449  fD0ZResPUpgSA ->SetName("D0ZResPUpgSA" );
450  }
451  if(resfileUpg->Get("D0ZResKSA" )) {
452  fD0ZResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResKSA" )));
453  fD0ZResKUpgSA ->SetName("D0ZResKUpgSA" );
454  }
455  if(resfileUpg->Get("D0ZResPiSA" )) {
456  fD0ZResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPiSA" )));
457  fD0ZResPiUpgSA ->SetName("D0ZResPiUpgSA" );
458  }
459  if(resfileUpg->Get("D0ZResESA" )) {
460  fD0ZResEUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResESA" )));
461  fD0ZResEUpgSA ->SetName("D0ZResEUpgSA" );
462  }
463  if(resfileUpg->Get("Pt1ResPSA" )) {
464  fPt1ResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPSA" )));
465  fPt1ResPUpgSA ->SetName("Pt1ResPUpgSA" );
466  }
467  if(resfileUpg->Get("Pt1ResKSA" )) {
468  fPt1ResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResKSA" )));
469  fPt1ResKUpgSA ->SetName("Pt1ResKUpgSA" );
470  }
471  if(resfileUpg->Get("Pt1ResPiSA" )) {
472  fPt1ResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPiSA" )));
473  fPt1ResPiUpgSA ->SetName("Pt1ResPiUpgSA" );
474  }
475  if(resfileUpg->Get("Pt1ResESA" )) {
476  fPt1ResEUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResESA" )));
477  fPt1ResEUpgSA ->SetName("Pt1ResEUpgSA" );
478  }
479  delete resfileUpg;
480  }
481 
482  DefineOutput(1,TNtuple::Class());
483 }
484 
486  //
487  // Destructor.
488  //
489  delete fDebugOutput;
490 }
491 
493  //
494  // Creation of user output objects.
495  //
496  fDebugOutput=new TList();
497  fDebugOutput->SetOwner();
498  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");
499  fDebugVars=new Float_t[fDebugNtuple->GetNvar()];
500 
501  fDebugOutput->Add(fDebugNtuple );
502 
511  for(Int_t j=0; j<2; j++){
512  for(Int_t i=0; i<4; i++){
513  if(fD0RPMeanPCur[j][i]) fDebugOutput->Add(fD0RPMeanPCur[j][i]);
514  if(fD0RPMeanKCur[j][i]) fDebugOutput->Add(fD0RPMeanKCur[j][i]);
515  if(fD0RPMeanPiCur[j][i]) fDebugOutput->Add(fD0RPMeanPiCur[j][i]);
516  if(fD0RPMeanECur[j][i]) fDebugOutput->Add(fD0RPMeanECur[j][i]);
517  if(fD0RPMeanPUpg[j][i]) fDebugOutput->Add(fD0RPMeanPUpg[j][i]);
518  if(fD0RPMeanKUpg[j][i]) fDebugOutput->Add(fD0RPMeanKUpg[j][i]);
519  if(fD0RPMeanPiUpg[j][i]) fDebugOutput->Add(fD0RPMeanPiUpg[j][i]);
520  if(fD0RPMeanEUpg[j][i]) fDebugOutput->Add(fD0RPMeanEUpg[j][i]);
521  }
522  }
543 
544  PostData(1,fDebugOutput);
545 }
546 
548  //
549  // The event loop
550  //
551  AliAODEvent *ev=0x0;
552  AliESDEvent *evesd=0x0;
553  Double_t bz=0.;
554 
555  if(fIsAOD) {
556  if(!fRunInVertexing) {
557  ev=dynamic_cast<AliAODEvent*>(InputEvent());
558  } else {
559  if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());
560  }
561  if(!ev) return;
562  bz=ev->GetMagneticField();
563  }
564  else {
565  evesd = dynamic_cast<AliESDEvent*>(InputEvent());
566  if (!evesd) {
567  AliError("event not found. Nothing done!");
568  return;
569  }
570  bz=evesd->GetMagneticField();
571  }
572 
573  if(fIsAOD) {
574 
575  // first loop on candidates to fill them in case of reduced AODs
576  // this is done to have the same behaviour of the improver with full (pp, p-Pb) and recuced (Pb-Pb) candidates
578 
579  // D0->Kpi
580  TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
581  if (array2Prong) {
582  for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
583  AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
584  vHF->FillRecoCand(ev,(AliAODRecoDecayHF2Prong*)decay);
585  }
586  }
587  // Dstar->Kpipi
588  TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
589  if (arrayCascade) {
590  for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
591  AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
592  vHF->FillRecoCasc(ev,((AliAODRecoCascadeHF*)decayDstar),kTRUE);
593  }
594  }
595  // Three prong
596  TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
597  if (array3Prong) {
598  for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
599  AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
600  vHF->FillRecoCand(ev,(AliAODRecoDecayHF3Prong*)decay);
601  }
602  }
603 
604 
605  // Smear all tracks
606  fMCs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
607  if (!fMCs) return;
608  if (fImproveTracks) {
609  for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
610  AliAODTrack * trk = dynamic_cast<AliAODTrack*>(ev->GetTrack(itrack));
611  if(!trk) AliFatal("Not a standard AOD");
612  SmearTrack(trk,bz);
613  }
614  }
615 
616  // TODO: recalculated primary vertex
617  AliVVertex *primaryVertex=ev->GetPrimaryVertex();
618 
619 
620  // Recalculate all candidates
621  // D0->Kpi
622  if (array2Prong) {
623  for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
624  AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
625  if(!vHF->FillRecoCand(ev,(AliAODRecoDecayHF2Prong*)decay))continue;
626 
627  // recalculate vertices
628  AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
629 
630 
631  AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
632  AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
633 
634  TObjArray ta12;
635 
636  ta12.Add(&et1); ta12.Add(&et2);
637  AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
638 
639 
640  // update secondary vertex
641  Double_t pos[3];
642  Double_t covpos[6];
643  v12->GetXYZ(pos);
644  v12->GetCovMatrix(covpos);
645  decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
646  if(fUpdateSecVertCovMat) decay->GetSecondaryVtx()->SetCovMatrix(covpos);
647  decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
648 
649  // update d0
650  Double_t d0z0[2],covd0z0[3];
651  Double_t d0[2],d0err[2];
652  et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
653  d0[0]=d0z0[0];
654  d0err[0] = TMath::Sqrt(covd0z0[0]);
655  et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
656  d0[1]=d0z0[0];
657  d0err[1] = TMath::Sqrt(covd0z0[0]);
658  decay->Setd0Prongs(2,d0);
659  decay->Setd0errProngs(2,d0err);
660  //
661 
662 
663  Double_t xdummy=0.,ydummy=0.;
664  Double_t dca;
665  dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
666  decay->SetDCA(dca);
667 
668 
669 
670  Double_t px[2],py[2],pz[2];
671  for (Int_t i=0;i<2;++i) {
672  AliExternalTrackParam et;
673  et.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(i)));
674  et.PropagateToDCA(v12,bz,100.,d0z0,covd0z0);
675  px[i]=et.Px();
676  py[i]=et.Py();
677  pz[i]=et.Pz();
678  }
679  decay->SetPxPyPzProngs(2,px,py,pz);
680  delete v12;
681  }
682  }
683 
684 
685  // Dstar->Kpipi
686  if (arrayCascade) {
687  for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
688  AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
689  if(!vHF->FillRecoCasc(ev,((AliAODRecoCascadeHF*)decayDstar),kTRUE))continue;
690  //Get D0 from D*
692 
693  // recalculate vertices
694  //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
695 
696  //soft pion
697  AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
698 
699  //track D0
700  AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
701 
703 
704  // update d0
705  Double_t d0z0[2],covd0z0[3];
706  Double_t d01[2],d01err[2];
707 
708  //the D*
709  et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
710  d01[0]=d0z0[0];
711  d01err[0] = TMath::Sqrt(covd0z0[0]);
712  trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
713  d01[1]=d0z0[0];
714  d01err[1] = TMath::Sqrt(covd0z0[0]);
715  decayDstar->Setd0Prongs(2,d01);
716  decayDstar->Setd0errProngs(2,d01err);
717 
718  // delete v12;
719  delete trackD0; trackD0=NULL;
720 
721  // a run for D*
722  Double_t px1[2],py1[2],pz1[2];
723  for (Int_t i=0;i<2;++i) {
724  const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));
725  px1[i]=t1->Px();
726  py1[i]=t1->Py();
727  pz1[i]=t1->Pz();
728  }
729  decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
730 
731  }
732  }
733 
734 
735  // Three prong
736  if (array3Prong) {
737  for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
738  AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
739  if(!vHF->FillRecoCand(ev,(AliAODRecoDecayHF3Prong*)decay))continue;
740 
741  // recalculate vertices
742  AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
743  AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
744  AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
745  AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
746  TObjArray ta123,ta12,ta23;
747  ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
748  ta12. Add(&et1);ta12 .Add(&et2);
749  ta23 .Add(&et2);ta23 .Add(&et3);
750  AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
751  AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
752  AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
753 
754  // update secondary vertex
755  Double_t pos[3];
756  Double_t covpos[6];
757  v123->GetXYZ(pos);
758  v123->GetCovMatrix(covpos);
759  decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
760  if(fUpdateSecVertCovMat) decay->GetSecondaryVtx()->SetCovMatrix(covpos);
761  decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
762 
763  // update d0 for all progs
764  Double_t d0z0[2],covd0z0[3];
765  Double_t d0[3],d0err[3];
766  et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
767  d0[0]=d0z0[0];
768  d0err[0] = TMath::Sqrt(covd0z0[0]);
769  et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
770  d0[1]=d0z0[0];
771  d0err[1] = TMath::Sqrt(covd0z0[0]);
772  et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
773  d0[2]=d0z0[0];
774  d0err[2] = TMath::Sqrt(covd0z0[0]);
775  decay->Setd0Prongs (3,d0 );
776  decay->Setd0errProngs(3,d0err);
777  // TODO: setter missing
778 
779  // update dca for prong combinations
780  Double_t xdummy=0.,ydummy=0.;
781  Double_t dca[3];
782  dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
783  dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
784  dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
785  decay->SetDCAs(3,dca);
786 
787  // update sigmavertex = dispersion
788  Float_t sigmaV=v123->GetDispersion();
789  decay->SetSigmaVert(sigmaV);
790  // update dist12 and dist23
791  primaryVertex->GetXYZ(pos);
792  decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
793  +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
794  +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
795  decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
796  +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
797  +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
798 
799 
800  Double_t px[3],py[3],pz[3];
801  for (Int_t i=0;i<3;++i) {
802  AliExternalTrackParam et;
803  et.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(i)));
804  et.PropagateToDCA(v123,bz,100.,d0z0,covd0z0);
805  px[i]=et.Px();
806  py[i]=et.Py();
807  pz[i]=et.Pz();
808  }
809  decay->SetPxPyPzProngs(3,px,py,pz);
810 
811  delete v123;delete v12;delete v23;
812  }
813  }
814  delete vHF;
815 
816  } // end AOD
817  else {
818  //
819  // In case of ESD: only smear all tracks
820  //
821  if (!fMCEvent) return;
822  if (fImproveTracks) {
823  for(Int_t itrack=0;itrack<evesd->GetNumberOfTracks();++itrack) {
824  AliESDtrack * trk = dynamic_cast<AliESDtrack*>(evesd->GetTrack(itrack));
825  if(!trk) AliFatal("No a standard ESD");
826  SmearTrack(trk,bz);
827  }
828  }
829  }// end ESD
830 
831 }
832 
834  // Early exit, if this track has nothing in common with the ITS
835 
836  if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
837  return;
838 
839  // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))
840  if (TESTBIT(track->GetITSClusterMap(),7)) return;
841  //
842 
843 
844  // Get reconstructed track parameters
845  AliExternalTrackParam et; et.CopyFromVTrack(track);
846  Double_t *param=const_cast<Double_t*>(et.GetParameter());
847 
848 
849  Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
850 
851  // Get MC info
852  Int_t imc=track->GetLabel();
853  if (imc<=0) return;
854  Double_t mcx[3];
855  Double_t mcp[3];
856  Double_t mccv[36]={0.};
857  Short_t mcc;
858  const AliVParticle *mc= 0x0;
859  if(fIsAOD) {
860  mc = dynamic_cast<AliVParticle*>(fMCs->At(imc));
861  }
862  else {
863  mc = dynamic_cast<AliVParticle*>(fMCEvent->GetTrack(imc));
864  }
865  if(!mc) return;
866  mc->XvYvZv(mcx);
867  mc->PxPyPz(mcp);
868  mcc=mc->Charge();
869  AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
870  const Double_t *parammc=mct.GetParameter();
871 //TODO: const Double_t *covermc=mct.GetCovariance();
872  AliVertex vtx(mcx,1.,1);
873 
874  // Correct reference points and frames according to MC
875  // TODO: B-Field correct?
876  // TODO: failing propagation....
877  et.PropagateToDCA(&vtx,track->GetBz(),10.);
878  et.Rotate(mct.GetAlpha());
879 
880  // Select appropriate smearing functions
881  Double_t ptmc=TMath::Abs(mc->Pt());
882  Double_t phimc=mc->Phi();
883  Int_t phiBin=PhiBin(phimc);
884  Int_t magfield=0;
885  if(bz<0.) magfield=0;
886  else if(bz>0.)magfield=1;
887  Double_t sd0rpn=0.;
888  Double_t sd0mrpn=0.;
889  Double_t sd0zn =0.;
890  Double_t spt1n =0.;
891  Double_t sd0rpo=0.;
892  Double_t sd0mrpo=0.;
893  Double_t sd0zo =0.;
894  Double_t spt1o =0.;
895  Double_t pullcorr=1.;
896  Int_t pdgcode = 0;
897  if(fIsAOD){
898  const AliAODMCParticle *amcpart = static_cast<const AliAODMCParticle *>(mc);
899  if(!amcpart) return;
900  pdgcode = amcpart->GetPdgCode();
901  } else {
902  const AliMCParticle *emcpart = static_cast<const AliMCParticle *>(mc);
903  if(!emcpart) return;
904  pdgcode = emcpart->PdgCode();
905  }
906 
907  switch (pdgcode) {
908  case 2212: case -2212:
910  sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);
911  spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);
913  sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);
914  spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);
915  sd0mrpo=EvalGraph(ptmc,fD0RPMeanPCur[magfield][phiBin],fD0RPMeanPCur[magfield][phiBin]);
916  sd0mrpn=EvalGraph(ptmc,fD0RPMeanPUpg[magfield][phiBin],fD0RPMeanPUpg[magfield][phiBin]);
918  break;
919  case 321: case -321:
921  sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);
922  spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);
924  sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);
925  spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);
926  sd0mrpo=EvalGraph(ptmc,fD0RPMeanKCur[magfield][phiBin],fD0RPMeanKCur[magfield][phiBin]);
927  sd0mrpn=EvalGraph(ptmc,fD0RPMeanKUpg[magfield][phiBin],fD0RPMeanKUpg[magfield][phiBin]);
929  break;
930  case 211: case -211:
937  sd0mrpo=EvalGraph(ptmc,fD0RPMeanPiCur[magfield][phiBin],fD0RPMeanPiCur[magfield][phiBin]);
938  sd0mrpn=EvalGraph(ptmc,fD0RPMeanPiUpg[magfield][phiBin],fD0RPMeanPiUpg[magfield][phiBin]);
940  break;
941  case 11: case -11:
943  sd0zo =EvalGraph(ptmc,fD0ZResECur,fD0ZResECurSA);
944  spt1o =EvalGraph(ptmc,fPt1ResECur,fPt1ResECurSA);
946  sd0zn =EvalGraph(ptmc,fD0ZResEUpg,fD0ZResEUpgSA);
947  spt1n =EvalGraph(ptmc,fPt1ResEUpg,fPt1ResEUpgSA);
948  sd0mrpo=EvalGraph(ptmc,fD0RPMeanECur[magfield][phiBin],fD0RPMeanECur[magfield][phiBin]);
949  sd0mrpn=EvalGraph(ptmc,fD0RPMeanEUpg[magfield][phiBin],fD0RPMeanEUpg[magfield][phiBin]);
951  break;
952  default:
953  return;
954  }
955 
956  // Use the same units (i.e. cm and GeV/c)! TODO: pt!
957  sd0rpo*=1.e-4;
958  sd0zo *=1.e-4;
959  sd0rpn*=1.e-4;
960  sd0zn *=1.e-4;
961  sd0mrpo*=1.e-4;
962  sd0mrpn*=1.e-4;
963 
964  // Apply the smearing
965  Double_t d0zo =param [1];
966  Double_t d0zmc =parammc[1];
967  Double_t d0rpo =param [0];
968  Double_t d0rpmc=parammc[0];
969  Double_t pt1o =param [4];
970  Double_t pt1mc =parammc[4];
971  Double_t dd0zo =d0zo-d0zmc;
972  Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
973  Double_t d0zn =d0zmc+dd0zn;
974  Double_t dd0rpo=d0rpo-d0rpmc;
975  Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
976  Double_t dd0mrpn=TMath::Abs(sd0mrpn)-TMath::Abs(sd0mrpo);
977  Double_t d0rpn =d0rpmc+dd0rpn-dd0mrpn;
978  Double_t d0zoinsigma = 0.;
979  if(covar[0] > 0.) d0zoinsigma = d0zo/TMath::Sqrt(covar[2]);
980  Double_t d0rpoinsigma = 0.;
981  if(covar[2] > 0.) d0rpoinsigma = d0rpo/TMath::Sqrt(covar[0]);
982 
983  if(fMimicData){
984  dd0mrpn=sd0mrpn-sd0mrpo;
985  d0rpn =d0rpmc+dd0rpn+dd0mrpn;
986  }
987 
988  Double_t dpt1o =pt1o-pt1mc;
989  Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
990  Double_t pt1n =pt1mc+dpt1n;
991  param[0]=d0rpn;
992  param[1]=d0zn ;
993  param[4]=pt1n ;
994 
995  //cov matrix update
996  if(fUpdateSTCovMatrix){
997  if(sd0rpo>0.) covar[0]*=(sd0rpn/sd0rpo)*(sd0rpn/sd0rpo);//yy
998  if(sd0zo>0. && sd0rpo>0.)covar[1]*=(sd0rpn/sd0rpo)*(sd0zn/sd0zo);//yz
999  if(sd0zo>0.) covar[2]*=(sd0zn/sd0zo)*(sd0zn/sd0zo);//zz
1000  if(sd0rpo>0.) covar[3]*=(sd0rpn/sd0rpo);//yl
1001  if(sd0zo>0.) covar[4]*=(sd0zn/sd0zo);//zl
1002  if(sd0rpo>0.) covar[6]*=(sd0rpn/sd0rpo);//ysenT
1003  if(sd0zo>0.) covar[7]*=(sd0zn/sd0zo);//zsenT
1004  if(sd0rpo>0. && spt1o>0.)covar[10]*=(sd0rpn/sd0rpo)*(spt1n/spt1o);//ypt
1005  if(sd0zo>0. && spt1o>0.) covar[11]*=(sd0zn/sd0zo)*(spt1n/spt1o);//zpt
1006  if(spt1o>0.) covar[12]*=(spt1n/spt1o);//sinPhipt
1007  if(spt1o>0.) covar[13]*=(spt1n/spt1o);//tanTpt
1008  if(spt1o>0.) covar[14]*=(spt1n/spt1o)*(spt1n/spt1o);//ptpt
1009  }
1010  if(fUpdatePulls){
1011 
1012  covar[0]*=pullcorr*pullcorr;//yy
1013  covar[1]*=pullcorr;//yz
1014  covar[3]*=pullcorr;//yl
1015  covar[6]*=pullcorr;//ysenT
1016  covar[10]*=pullcorr;//ypt
1017 
1018  }
1019 
1020  Double_t d0zninsigma = 0.;
1021  if(covar[0] > 0.) d0zninsigma = d0zn/TMath::Sqrt(covar[2]);
1022  Double_t d0rpninsigma = 0.;
1023  if(covar[2] > 0.) d0rpninsigma = d0rpn/TMath::Sqrt(covar[0]);
1024 
1025  // Copy the smeared parameters to the AOD track
1026  Double_t x[3];
1027  Double_t p[3];
1028  et.GetXYZ(x);
1029  et.GetPxPyPz(p);
1030  Double_t cv[21];
1031  et.GetCovarianceXYZPxPyPz(cv);
1032  if(fIsAOD) {
1033  AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
1034  aodtrack->SetPosition(x,kFALSE);
1035  aodtrack->SetP(p,kTRUE);
1036  aodtrack->SetCovMatrix(cv);
1037  // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
1038  UChar_t itsClusterMap = aodtrack->GetITSClusterMap();
1039  SETBIT(itsClusterMap,7);
1040  aodtrack->SetITSClusterMap(itsClusterMap);
1041  //
1042 
1043  } else {
1044  AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
1045  Short_t sign = esdtrack->Charge();
1046  esdtrack->Set(x,p,cv,sign);
1047  esdtrack->RelateToVVertex(InputEvent()->GetPrimaryVertex(), bz,100.);
1048  // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
1049  UChar_t itsClusterMap = esdtrack->GetITSClusterMap();
1050  SETBIT(itsClusterMap,7);
1051  esdtrack->SetITSClusterMap(itsClusterMap);
1052  }
1053 
1054 
1055  // write out debug infos
1056  if (fDebugNtuple->GetEntries()<fNDebug) {
1057  Int_t idbg=0;
1058  fDebugVars[idbg++]=pdgcode;
1059  fDebugVars[idbg++]=ptmc ;
1060  fDebugVars[idbg++]=d0rpo ;
1061  fDebugVars[idbg++]=d0zo ;
1062  fDebugVars[idbg++]=pt1o ;
1063  fDebugVars[idbg++]=sd0rpo;
1064  fDebugVars[idbg++]=sd0zo ;
1065  fDebugVars[idbg++]=spt1o ;
1066  fDebugVars[idbg++]=d0rpn ;
1067  fDebugVars[idbg++]=d0zn ;
1068  fDebugVars[idbg++]=pt1n ;
1069  fDebugVars[idbg++]=sd0rpn;
1070  fDebugVars[idbg++]=sd0zn ;
1071  fDebugVars[idbg++]=spt1n ;
1072  fDebugVars[idbg++]=d0rpmc;
1073  fDebugVars[idbg++]=d0zmc ;
1074  fDebugVars[idbg++]=pt1mc ;
1075  fDebugVars[idbg++]=pullcorr ;
1076  fDebugVars[idbg++]=d0zoinsigma ;
1077  fDebugVars[idbg++]=d0zninsigma ;
1078  fDebugVars[idbg++]=d0rpoinsigma ;
1079  fDebugVars[idbg++]=d0rpninsigma ;
1080  fDebugNtuple->Fill(fDebugVars);
1081  PostData(1,fDebugOutput);
1082  }
1083 }
1084 
1085 AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
1086  //
1087  // Helper function to recalculate a vertex.
1088  //
1089 
1090  static UShort_t ids[]={1,2,3}; //TODO: unsave...
1091  AliVertexerTracks vertexer(bField);
1092  vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
1093  AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
1094  return vertex;
1095 }
1096 
1098  //
1099  // Evaluates a TGraph without linear extrapolation. Instead the last
1100  // valid point of the graph is used when out of range.
1101  // The function assumes an ascending order of X.
1102  //
1103 
1104  if(!graph) return 0.;
1105 
1106  // TODO: find a pretty solution for this:
1107  Int_t n =graph->GetN();
1108  Double_t xmin=graph->GetX()[0 ];
1109  Double_t xmax=graph->GetX()[n-1];
1110  if (x<xmin) {
1111  if(!graphSA) return graph->Eval(xmin);
1112  Double_t xminSA=graphSA->GetX()[0];
1113  if(x<xminSA) return graphSA->Eval(xminSA);
1114  return graphSA->Eval(x);
1115  }
1116  if (x>xmax) return graph->Eval(xmax);
1117  return graph->Eval(x);
1118 }
1119 
1120 //________________________________________________________________________
1121 
1123  Double_t pi=TMath::Pi();
1124  if(phi>2.*pi || phi<0.) return -1;
1125  if((phi<=(pi/4.)) || (phi>7.*(pi/4.))) return 0;
1126  if((phi>(pi/4.)) && (phi<=3.*(pi/4.))) return 1;
1127  if((phi>3.*(pi/4.)) && (phi<=5.*(pi/4.))) return 2;
1128  if((phi>(5.*pi/4.)) && (phi<=7.*(pi/4.))) return 3;
1129  return -1;
1130 }
1131 
1132 ClassImp(AliAnalysisTaskSEImproveITS);
1133 
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