AliPhysics  8bb951a (8bb951a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 "AliVertexerTracks.h"
27 #include "AliAODEvent.h"
28 #include "AliAODTrack.h"
29 #include "AliAODMCParticle.h"
30 #include "AliExternalTrackParam.h"
33 #include "AliAODRecoCascadeHF.h"
34 #include "AliNeutralTrackParam.h"
36 
37 //
38 // Implementation of the "hybrid-approach" for ITS upgrade studies.
39 // The tastk smears the track parameters according to estimations
40 // from single-track upgrade studies. Afterwards it recalculates
41 // the parameters of the reconstructed decays.
42 //
43 // WARNING: This will affect all tasks in a train after this one
44 // (which is typically desired, though).
45 //
46 
48  :AliAnalysisTaskSE(),
49  fD0ZResPCur (0),
50  fD0ZResKCur (0),
51  fD0ZResPiCur (0),
52  fD0RPResPCur (0),
53  fD0RPResKCur (0),
54  fD0RPResPiCur(0),
55  fPt1ResPCur (0),
56  fPt1ResKCur (0),
57  fPt1ResPiCur (0),
58  fD0ZResPUpg (0),
59  fD0ZResKUpg (0),
60  fD0ZResPiUpg (0),
61  fD0RPResPUpg (0),
62  fD0RPResKUpg (0),
63  fD0RPResPiUpg(0),
64  fPt1ResPUpg (0),
65  fPt1ResKUpg (0),
66  fPt1ResPiUpg (0),
67  fD0ZResPCurSA (0),
68  fD0ZResKCurSA (0),
69  fD0ZResPiCurSA (0),
70  fD0RPResPCurSA (0),
71  fD0RPResKCurSA (0),
72  fD0RPResPiCurSA(0),
73  fPt1ResPCurSA (0),
74  fPt1ResKCurSA (0),
75  fPt1ResPiCurSA (0),
76  fD0ZResPUpgSA (0),
77  fD0ZResKUpgSA (0),
78  fD0ZResPiUpgSA (0),
79  fD0RPResPUpgSA (0),
80  fD0RPResKUpgSA (0),
81  fD0RPResPiUpgSA(0),
82  fPt1ResPUpgSA (0),
83  fPt1ResKUpgSA (0),
84  fPt1ResPiUpgSA (0),
85  fRunInVertexing(kFALSE),
86  fImproveTracks(kTRUE),
87  fDebugOutput (0),
88  fDebugNtuple (0),
89  fDebugVars (0),
90  fNDebug (0)
91 {
92  //
93  // Default constructor.
94  //
95 }
96 
98  const char *resfileCurURI,
99  const char *resfileUpgURI,
100  Bool_t isRunInVertexing,
101  Int_t ndebug)
102  :AliAnalysisTaskSE(name),
103  fD0ZResPCur (0),
104  fD0ZResKCur (0),
105  fD0ZResPiCur (0),
106  fD0RPResPCur (0),
107  fD0RPResKCur (0),
108  fD0RPResPiCur(0),
109  fPt1ResPCur (0),
110  fPt1ResKCur (0),
111  fPt1ResPiCur (0),
112  fD0ZResPUpg (0),
113  fD0ZResKUpg (0),
114  fD0ZResPiUpg (0),
115  fD0RPResPUpg (0),
116  fD0RPResKUpg (0),
117  fD0RPResPiUpg(0),
118  fPt1ResPUpg (0),
119  fPt1ResKUpg (0),
120  fPt1ResPiUpg (0),
121  fD0ZResPCurSA (0),
122  fD0ZResKCurSA (0),
123  fD0ZResPiCurSA (0),
124  fD0RPResPCurSA (0),
125  fD0RPResKCurSA (0),
126  fD0RPResPiCurSA(0),
127  fPt1ResPCurSA (0),
128  fPt1ResKCurSA (0),
129  fPt1ResPiCurSA (0),
130  fD0ZResPUpgSA (0),
131  fD0ZResKUpgSA (0),
132  fD0ZResPiUpgSA (0),
133  fD0RPResPUpgSA (0),
134  fD0RPResKUpgSA (0),
135  fD0RPResPiUpgSA(0),
136  fPt1ResPUpgSA (0),
137  fPt1ResKUpgSA (0),
138  fPt1ResPiUpgSA (0),
139  fRunInVertexing(isRunInVertexing),
140  fImproveTracks(kTRUE),
141  fDebugOutput (0),
142  fDebugNtuple (0),
143  fDebugVars (0),
144  fNDebug (ndebug)
145 {
146  //
147  // Constructor to be used to create the task.
148  // The the URIs specify the resolution files to be used.
149  // They are expected to contain TGraphs with the resolutions
150  // for the current and the upgraded ITS (see code for details).
151  // One may also specify for how many tracks debug information
152  // is written to the output.
153  //
154  TFile *resfileCur=TFile::Open(resfileCurURI);
155  fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
156  fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
157  fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
158  fD0ZResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP" )));
159  fD0ZResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK" )));
160  fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
161  fPt1ResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP" )));
162  fPt1ResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK" )));
163  fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
164  fD0RPResPCur ->SetName("D0RPResPCur" );
165  fD0RPResKCur ->SetName("D0RPResKCur" );
166  fD0RPResPiCur->SetName("D0RPResPiCur");
167  fD0ZResPCur ->SetName("D0ZResPCur" );
168  fD0ZResKCur ->SetName("D0ZResKCur" );
169  fD0ZResPiCur ->SetName("D0ZResPiCur" );
170  fPt1ResPCur ->SetName("Pt1ResPCur" );
171  fPt1ResKCur ->SetName("Pt1ResKCur" );
172  fPt1ResPiCur ->SetName("Pt1ResPiCur" );
173  fD0RPResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPSA" )));
174  fD0RPResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResKSA" )));
175  fD0RPResPiCurSA=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPiSA")));
176  fD0ZResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPSA" )));
177  fD0ZResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResKSA" )));
178  fD0ZResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPiSA" )));
179  fPt1ResPCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPSA" )));
180  fPt1ResKCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResKSA" )));
181  fPt1ResPiCurSA =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPiSA" )));
182  fD0RPResPCurSA ->SetName("D0RPResPCurSA" );
183  fD0RPResKCurSA ->SetName("D0RPResKCurSA" );
184  fD0RPResPiCurSA->SetName("D0RPResPiCurSA");
185  fD0ZResPCurSA ->SetName("D0ZResPCurSA" );
186  fD0ZResKCurSA ->SetName("D0ZResKCurSA" );
187  fD0ZResPiCurSA ->SetName("D0ZResPiCurSA" );
188  fPt1ResPCurSA ->SetName("Pt1ResPCurSA" );
189  fPt1ResKCurSA ->SetName("Pt1ResKCurSA" );
190  fPt1ResPiCurSA ->SetName("Pt1ResPiCurSA" );
191  delete resfileCur;
192  TFile *resfileUpg=TFile::Open(resfileUpgURI );
193  fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
194  fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
195  fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
196  fD0ZResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP" )));
197  fD0ZResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK" )));
198  fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
199  fPt1ResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP" )));
200  fPt1ResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK" )));
201  fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
202  fD0RPResPUpg ->SetName("D0RPResPUpg" );
203  fD0RPResKUpg ->SetName("D0RPResKUpg" );
204  fD0RPResPiUpg->SetName("D0RPResPiUpg");
205  fD0ZResPUpg ->SetName("D0ZResPUpg" );
206  fD0ZResKUpg ->SetName("D0ZResKUpg" );
207  fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
208  fPt1ResPUpg ->SetName("Pt1ResPUpg" );
209  fPt1ResKUpg ->SetName("Pt1ResKUpg" );
210  fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
211  fD0RPResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPSA" )));
212  fD0RPResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResKSA" )));
213  fD0RPResPiUpgSA=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPiSA")));
214  fD0ZResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPSA" )));
215  fD0ZResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResKSA" )));
216  fD0ZResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPiSA" )));
217  fPt1ResPUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPSA" )));
218  fPt1ResKUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResKSA" )));
219  fPt1ResPiUpgSA =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPiSA" )));
220  fD0RPResPUpgSA ->SetName("D0RPResPUpgSA" );
221  fD0RPResKUpgSA ->SetName("D0RPResKUpgSA" );
222  fD0RPResPiUpgSA->SetName("D0RPResPiUpgSA");
223  fD0ZResPUpgSA ->SetName("D0ZResPUpgSA" );
224  fD0ZResKUpgSA ->SetName("D0ZResKUpgSA" );
225  fD0ZResPiUpgSA ->SetName("D0ZResPiUpgSA" );
226  fPt1ResPUpgSA ->SetName("Pt1ResPUpgSA" );
227  fPt1ResKUpgSA ->SetName("Pt1ResKUpgSA" );
228  fPt1ResPiUpgSA ->SetName("Pt1ResPiUpgSA" );
229  delete resfileUpg;
230 
231  DefineOutput(1,TNtuple::Class());
232 }
233 
235  //
236  // Destructor.
237  //
238  delete fDebugOutput;
239 }
240 
242  //
243  // Creation of user output objects.
244  //
245  fDebugOutput=new TList();
246  fDebugOutput->SetOwner();
247  fDebugNtuple=new TNtuple("fDebugNtuple","Smearing","pdg:ptmc:d0rpo:d0zo:pt1o:sd0rpo:sd0zo:spt1o:d0rpn:d0zn:pt1n:sd0rpn:sd0zn:spt1n:d0rpmc:d0zmc:pt1mc");
248  fDebugVars=new Float_t[fDebugNtuple->GetNvar()];
249 
250  fDebugOutput->Add(fDebugNtuple );
251 
252  fDebugOutput->Add(fD0RPResPCur );
253  fDebugOutput->Add(fD0RPResKCur );
255  fDebugOutput->Add(fD0ZResPCur );
256  fDebugOutput->Add(fD0ZResKCur );
257  fDebugOutput->Add(fD0ZResPiCur );
258  fDebugOutput->Add(fPt1ResPCur );
259  fDebugOutput->Add(fPt1ResKCur );
260  fDebugOutput->Add(fPt1ResPiCur );
261  fDebugOutput->Add(fD0RPResPUpg );
262  fDebugOutput->Add(fD0RPResKUpg );
264  fDebugOutput->Add(fD0ZResPUpg );
265  fDebugOutput->Add(fD0ZResKUpg );
266  fDebugOutput->Add(fD0ZResPiUpg );
267  fDebugOutput->Add(fPt1ResPUpg );
268  fDebugOutput->Add(fPt1ResKUpg );
269  fDebugOutput->Add(fPt1ResPiUpg );
270 
271  PostData(1,fDebugOutput);
272 }
273 
275  //
276  // The event loop
277  //
278  AliAODEvent *ev=0x0;
279  if(!fRunInVertexing) {
280  ev=dynamic_cast<AliAODEvent*>(InputEvent());
281  } else {
282  if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());
283  }
284  if(!ev) return;
285  Double_t bz=ev->GetMagneticField();
286 
287 
288 
289 
290  // Smear all tracks
291  TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
292  if (!mcs) return;
293  if (fImproveTracks) {
294  for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
295  AliAODTrack * trk = dynamic_cast<AliAODTrack*>(ev->GetTrack(itrack));
296  if(!trk) AliFatal("Not a standard AOD");
297  SmearTrack(trk,mcs);
298  }
299  }
300 
301  // TODO: recalculated primary vertex
302  AliVVertex *primaryVertex=ev->GetPrimaryVertex();
303 
304  // Recalculate all candidates
305  // D0->Kpi
306  TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
307  if (array2Prong) {
308  for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
309  AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
310 
311  // recalculate vertices
312  AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
313 
314 
315  AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
316  AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
317 
318  TObjArray ta12;
319 
320  ta12.Add(&et1); ta12.Add(&et2);
321  AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
322 
323 
324  // update secondary vertex
325  Double_t pos[3];
326  v12->GetXYZ(pos);
327 
328  decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
329  decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
330 
332 
333  // update d0
334  Double_t d0z0[2],covd0z0[3];
335  Double_t d0[2],d0err[2];
336  et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
337  d0[0]=d0z0[0];
338  d0err[0] = TMath::Sqrt(covd0z0[0]);
339  et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
340  d0[1]=d0z0[0];
341  d0err[1] = TMath::Sqrt(covd0z0[0]);
342  decay->Setd0Prongs(2,d0);
343  decay->Setd0errProngs(2,d0err);
344  //
345 
346 
347  Double_t xdummy=0.,ydummy=0.;
348  Double_t dca;
349  dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
350  decay->SetDCA(dca);
351 
352 
353  delete v12;
354 
355  Double_t px[2],py[2],pz[2];
356  for (Int_t i=0;i<2;++i) {
357  const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
358  px[i]=t->Px();
359  py[i]=t->Py();
360  pz[i]=t->Pz();
361  }
362  decay->SetPxPyPzProngs(2,px,py,pz);
363  }
364  }
365 
366 
367  // Dstar->Kpipi
368  TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
369 
370  if (arrayCascade) {
371  for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
372  AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
373  //Get D0 from D*
375 
376  // recalculate vertices
377  //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
378 
379  //soft pion
380  AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
381 
382  //track D0
383  AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
384 
386 
387  // update d0
388  Double_t d0z0[2],covd0z0[3];
389  Double_t d01[2],d01err[2];
390 
391  //the D*
392  et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
393  d01[0]=d0z0[0];
394  d01err[0] = TMath::Sqrt(covd0z0[0]);
395  trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
396  d01[1]=d0z0[0];
397  d01err[1] = TMath::Sqrt(covd0z0[0]);
398  decayDstar->Setd0Prongs(2,d01);
399  decayDstar->Setd0errProngs(2,d01err);
400 
401  // delete v12;
402  delete trackD0; trackD0=NULL;
403 
404  // a run for D*
405  Double_t px1[2],py1[2],pz1[2];
406  for (Int_t i=0;i<2;++i) {
407  const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));
408  px1[i]=t1->Px();
409  py1[i]=t1->Py();
410  pz1[i]=t1->Pz();
411  }
412  decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
413 
414  }
415  }
416 
417 
418  // Three prong
419  TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
420  if (array3Prong) {
421  for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
422  AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
423 
424  // recalculate vertices
425  AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
426  AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
427  AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
428  AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
429  TObjArray ta123,ta12,ta23;
430  ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
431  ta12. Add(&et1);ta12 .Add(&et2);
432  ta23 .Add(&et2);ta23 .Add(&et3);
433  AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
434  AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
435  AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
436 
437  // update secondary vertex
438  Double_t pos[3];
439  v123->GetXYZ(pos);
440  decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
441  decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
442  //TODO: covariance matrix
443 
444  // update d0 for all progs
445  Double_t d0z0[2],covd0z0[3];
446  Double_t d0[3],d0err[3];
447  et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
448  d0[0]=d0z0[0];
449  d0err[0] = TMath::Sqrt(covd0z0[0]);
450  et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
451  d0[1]=d0z0[0];
452  d0err[1] = TMath::Sqrt(covd0z0[0]);
453  et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
454  d0[2]=d0z0[0];
455  d0err[2] = TMath::Sqrt(covd0z0[0]);
456  decay->Setd0Prongs (3,d0 );
457  decay->Setd0errProngs(3,d0err);
458  // TODO: setter missing
459 
460  // update dca for prong combinations
461  Double_t xdummy=0.,ydummy=0.;
462  Double_t dca[3];
463  dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
464  dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
465  dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
466  decay->SetDCAs(3,dca);
467 
468  // update dist12 and dist23
469  primaryVertex->GetXYZ(pos);
470  decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
471  +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
472  +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
473  decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
474  +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
475  +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
476 
477  delete v123;delete v12;delete v23;
478 
479  Double_t px[3],py[3],pz[3];
480  for (Int_t i=0;i<3;++i) {
481  const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
482  px[i]=t->Px();
483  py[i]=t->Py();
484  pz[i]=t->Pz();
485  }
486  decay->SetPxPyPzProngs(3,px,py,pz);
487  }
488  }
489 }
490 
491 void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArray *mcs) {
492  // Early exit, if this track has nothing in common with the ITS
493 
494  if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
495  return;
496 
497  // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))
498  if (TESTBIT(track->GetITSClusterMap(),7)) return;
499  //
500 
501 
502  // Get reconstructed track parameters
503  AliExternalTrackParam et; et.CopyFromVTrack(track);
504  Double_t *param=const_cast<Double_t*>(et.GetParameter());
505 //TODO: Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
506 
507  // Get MC info
508  Int_t imc=track->GetLabel();
509  if (imc<=0) return;
510  const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
511  Double_t mcx[3];
512  Double_t mcp[3];
513  Double_t mccv[36]={0.};
514  Short_t mcc;
515  mc->XvYvZv(mcx);
516  mc->PxPyPz(mcp);
517  mcc=mc->Charge();
518  AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
519  const Double_t *parammc=mct.GetParameter();
520 //TODO: const Double_t *covermc=mct.GetCovariance();
521  AliVertex vtx(mcx,1.,1);
522 
523  // Correct reference points and frames according to MC
524  // TODO: B-Field correct?
525  // TODO: failing propagation....
526  et.PropagateToDCA(&vtx,track->GetBz(),10.);
527  et.Rotate(mct.GetAlpha());
528 
529  // Select appropriate smearing functions
530  Double_t ptmc=TMath::Abs(mc->Pt());
531  Double_t sd0rpn=0.;
532  Double_t sd0zn =0.;
533  Double_t spt1n =0.;
534  Double_t sd0rpo=0.;
535  Double_t sd0zo =0.;
536  Double_t spt1o =0.;
537  switch (mc->GetPdgCode()) {
538  case 2212: case -2212:
540  sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);
541  spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);
543  sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);
544  spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);
545  break;
546  case 321: case -321:
548  sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);
549  spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);
551  sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);
552  spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);
553  break;
554  case 211: case -211:
561  break;
562  default:
563  return;
564  }
565 
566  // Use the same units (i.e. cm and GeV/c)! TODO: pt!
567  sd0rpo*=1.e-4;
568  sd0zo *=1.e-4;
569  sd0rpn*=1.e-4;
570  sd0zn *=1.e-4;
571 
572  // Apply the smearing
573  Double_t d0zo =param [1];
574  Double_t d0zmc =parammc[1];
575  Double_t d0rpo =param [0];
576  Double_t d0rpmc=parammc[0];
577  Double_t pt1o =param [4];
578  Double_t pt1mc =parammc[4];
579  Double_t dd0zo =d0zo-d0zmc;
580  Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
581  Double_t d0zn =d0zmc+dd0zn;
582  Double_t dd0rpo=d0rpo-d0rpmc;
583  Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
584  Double_t d0rpn =d0rpmc+dd0rpn;
585  Double_t dpt1o =pt1o-pt1mc;
586  Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
587  Double_t pt1n =pt1mc+dpt1n;
588  param[0]=d0rpn;
589  param[1]=d0zn ;
590  param[4]=pt1n ;
591 
592  // Copy the smeared parameters to the AOD track
593  Double_t x[3];
594  Double_t p[3];
595  et.GetXYZ(x);
596  et.GetPxPyPz(p);
597  track->SetPosition(x,kFALSE);
598  track->SetP(p,kTRUE);
599 
600 
601  // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
602  UChar_t itsClusterMap = track->GetITSClusterMap();
603  SETBIT(itsClusterMap,7);
604  track->SetITSClusterMap(itsClusterMap);
605  //
606 
607  // write out debug infos
608  if (fDebugNtuple->GetEntries()<fNDebug) {
609  Int_t idbg=0;
610  fDebugVars[idbg++]=mc->GetPdgCode();
611  fDebugVars[idbg++]=ptmc ;
612  fDebugVars[idbg++]=d0rpo ;
613  fDebugVars[idbg++]=d0zo ;
614  fDebugVars[idbg++]=pt1o ;
615  fDebugVars[idbg++]=sd0rpo;
616  fDebugVars[idbg++]=sd0zo ;
617  fDebugVars[idbg++]=spt1o ;
618  fDebugVars[idbg++]=d0rpn ;
619  fDebugVars[idbg++]=d0zn ;
620  fDebugVars[idbg++]=pt1n ;
621  fDebugVars[idbg++]=sd0rpn;
622  fDebugVars[idbg++]=sd0zn ;
623  fDebugVars[idbg++]=spt1n ;
624  fDebugVars[idbg++]=d0rpmc;
625  fDebugVars[idbg++]=d0zmc ;
626  fDebugVars[idbg++]=pt1mc ;
627  fDebugNtuple->Fill(fDebugVars);
628  PostData(1,fDebugOutput);
629  }
630 }
631 
632 AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
633  //
634  // Helper function to recalculate a vertex.
635  //
636 
637  static UShort_t ids[]={1,2,3}; //TODO: unsave...
638  AliVertexerTracks vertexer(bField);
639  vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
640  AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
641  return vertex;
642 }
643 
644 Double_t AliAnalysisTaskSEImproveITS::EvalGraph(Double_t x,const TGraph *graph,const TGraph *graphSA) const {
645  //
646  // Evaluates a TGraph without linear extrapolation. Instead the last
647  // valid point of the graph is used when out of range.
648  // The function assumes an ascending order of X.
649  //
650 
651  // TODO: find a pretty solution for this:
652  Int_t n =graph->GetN();
653  Double_t xmin=graph->GetX()[0 ];
654  Double_t xmax=graph->GetX()[n-1];
655  if (x<xmin) {
656  if(!graphSA) return graph->Eval(xmin);
657  Double_t xminSA=graphSA->GetX()[0];
658  if(x<xminSA) return graphSA->Eval(xminSA);
659  return graphSA->Eval(x);
660  }
661  if (x>xmax) return graph->Eval(xmax);
662  return graph->Eval(x);
663 }
664 
666 
TGraph * fPt1ResPiCur
old pt dep. 1/pt res. for kaons
TGraph * fD0ZResPUpg
old pt dep. 1/pt res. for pions
TGraph * fD0RPResKCur
old pt dep. d0 res. in rphi for protons
TGraph * fD0RPResPiUpgSA
new standalone pt dep. d0 res. in rphi for kaons
TGraph * fD0ZResKUpgSA
new standalone pt dep. d0 res. in z for protons
TGraph * fD0RPResPCurSA
old standalone pt dep. d0 res. in z for pions
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 pions
Double_t EvalGraph(Double_t x, const TGraph *graph, const TGraph *graphSA=0) const
Helper functions.
void SmearTrack(AliAODTrack *track, const TClonesArray *mcs)
TGraph * fD0ZResPiCurSA
old standalone pt dep. d0 res. in z for kaons
virtual void UserExec(Option_t *option)
TList * fDebugOutput
this is always kTRUE. kFALSE only if re-running on already improved AODs
TGraph * fPt1ResPiUpg
new pt dep. 1/pt res. for kaons
TGraph * fPt1ResPCurSA
old standalone pt dep. d0 res. in rphi for pions
TGraph * fPt1ResPiCurSA
old standalone pt dep. 1/pt res. for kaons
TGraph * fPt1ResPUpg
new pt dep. d0 res. in rphi for pions
TGraph * fD0RPResKCurSA
old standalone pt dep. d0 res. in rphi for protons
TGraph * fPt1ResKUpg
new pt dep. 1/pt res. for protons
TGraph * fPt1ResKCur
old pt dep. 1/pt res. for protons
TGraph * fD0RPResKUpgSA
new standalone pt dep. d0 res. in rphi for protons
TGraph * fPt1ResPUpgSA
new standalone pt dep. d0 res. in rphi for pions
AliAODTrack * GetBachelor() const
TGraph * fD0RPResPiCurSA
old standalone pt dep. d0 res. in rphi for kaons
TGraph * fD0ZResKCurSA
old standalone pt dep. d0 res. in z for protons
TGraph * fD0ZResPCurSA
new pt dep. 1/pt res. for pions
Float_t * fDebugVars
! variables to store as degug info
TGraph * fD0RPResPCur
old pt dep. d0 res. in z for pions
Bool_t fRunInVertexing
new standalone pt dep. 1/pt res. for pions
TGraph * fPt1ResPiUpgSA
new standalone pt dep. 1/pt res. for kaons
TGraph * fD0RPResPUpgSA
new standalone pt dep. d0 res. in z for pions
TGraph * fPt1ResKCurSA
old standalone pt dep. 1/pt res. for protons
Bool_t fImproveTracks
flag to run hybrid task before the vertexingHF task or in standard mode
TGraph * fD0ZResKCur
old pt dep. d0 res. in z for protons
void Setd0errProngs(Int_t nprongs, Double_t *d0)
TGraph * fPt1ResPCur
old pt dep. d0 res. in rphi for pions
TGraph * fD0RPResPUpg
new pt dep. d0 res. in z for pions
TGraph * fD0ZResPiUpg
new pt dep. d0 res. in z for kaons
TGraph * fD0ZResPiCur
old pt dep. d0 res. in z for kaons
AliAODRecoDecayHF2Prong * Get2Prong() const
ClassImp(AliAnalysisTaskSEImproveITS)
TGraph * fD0ZResPiUpgSA
new standalone pt dep. d0 res. in z for kaons
TGraph * fD0RPResPiUpg
new pt dep. d0 res. in rphi for kaons
AliESDVertex * RecalculateVertex(const AliVVertex *old, TObjArray *tracks, Double_t bField)
TNtuple * fDebugNtuple
! debug send on output slot 1
TGraph * fD0RPResKUpg
new pt dep. d0 res. in rphi for protons
TGraph * fD0ZResKUpg
new pt dep. d0 res. in z for protons