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