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