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