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