]> 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
end-of-line normalization
[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) {
295 SmearTrack(ev->GetTrack(itrack),mcs);
296 }
297 }
298
299 // TODO: recalculated primary vertex
300 AliVVertex *primaryVertex=ev->GetPrimaryVertex();
301
302 // Recalculate all candidates
303 // D0->Kpi
304 TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
305 if (array2Prong) {
306 for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
307 AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
308
309 // recalculate vertices
310 AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
311
312
313 AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
314 AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
315
316 TObjArray ta12;
317
318 ta12.Add(&et1); ta12.Add(&et2);
319 AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
320
321
322 // update secondary vertex
323 Double_t pos[3];
324 v12->GetXYZ(pos);
325
326 decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
327 decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
328
329 //!!!!TODO: covariance matrix
330
331 // update d0
332 Double_t d0z0[2],covd0z0[3];
333 Double_t d0[2],d0err[2];
334 et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
335 d0[0]=d0z0[0];
336 d0err[0] = TMath::Sqrt(covd0z0[0]);
337 et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
338 d0[1]=d0z0[0];
339 d0err[1] = TMath::Sqrt(covd0z0[0]);
340 decay->Setd0Prongs(2,d0);
341 decay->Setd0errProngs(2,d0err);
342 //
343
344
345 Double_t xdummy=0.,ydummy=0.;
346 Double_t dca;
347 dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
348 decay->SetDCA(dca);
349
350
351 delete v12;
352
353 Double_t px[2],py[2],pz[2];
354 for (Int_t i=0;i<2;++i) {
355 const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
356 px[i]=t->Px();
357 py[i]=t->Py();
358 pz[i]=t->Pz();
359 }
360 decay->SetPxPyPzProngs(2,px,py,pz);
361 }
362 }
363
364
365 // Dstar->Kpipi
366 TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
367
368 if (arrayCascade) {
369 for (Int_t icand=0;icand<arrayCascade->GetEntries();++icand) {
370 AliAODRecoCascadeHF *decayDstar=static_cast<AliAODRecoCascadeHF*>(arrayCascade->At(icand));
371 //Get D0 from D*
372 AliAODRecoDecayHF2Prong* decay=(AliAODRecoDecayHF2Prong*)decayDstar->Get2Prong();
373
374 // recalculate vertices
375 //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
376
377 //soft pion
378 AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
379
380 //track D0
381 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
382
383 //!!!!TODO: covariance matrix
384
385 // update d0
386 Double_t d0z0[2],covd0z0[3];
387 Double_t d01[2],d01err[2];
388
389 //the D*
390 et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
391 d01[0]=d0z0[0];
392 d01err[0] = TMath::Sqrt(covd0z0[0]);
393 trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
394 d01[1]=d0z0[0];
395 d01err[1] = TMath::Sqrt(covd0z0[0]);
396 decayDstar->Setd0Prongs(2,d01);
397 decayDstar->Setd0errProngs(2,d01err);
398
399 // delete v12;
400 delete trackD0; trackD0=NULL;
401
402 // a run for D*
403 Double_t px1[2],py1[2],pz1[2];
404 for (Int_t i=0;i<2;++i) {
405 const AliAODTrack *t1=static_cast<AliAODTrack*>(decayDstar->GetDaughter(i));
406 px1[i]=t1->Px();
407 py1[i]=t1->Py();
408 pz1[i]=t1->Pz();
409 }
410 decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
411
412 }
413 }
414
415
416 // Three prong
417 TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
418 if (array3Prong) {
419 for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
420 AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
421
422 // recalculate vertices
423 AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
424 AliExternalTrackParam et1; et1.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
425 AliExternalTrackParam et2; et2.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
426 AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
427 TObjArray ta123,ta12,ta23;
428 ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
429 ta12. Add(&et1);ta12 .Add(&et2);
430 ta23 .Add(&et2);ta23 .Add(&et3);
431 AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
432 AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
433 AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
434
435 // update secondary vertex
436 Double_t pos[3];
437 v123->GetXYZ(pos);
438 decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
439 decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
440 //TODO: covariance matrix
441
442 // update d0 for all progs
443 Double_t d0z0[2],covd0z0[3];
444 Double_t d0[3],d0err[3];
445 et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
446 d0[0]=d0z0[0];
447 d0err[0] = TMath::Sqrt(covd0z0[0]);
448 et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
449 d0[1]=d0z0[0];
450 d0err[1] = TMath::Sqrt(covd0z0[0]);
451 et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
452 d0[2]=d0z0[0];
453 d0err[2] = TMath::Sqrt(covd0z0[0]);
454 decay->Setd0Prongs (3,d0 );
455 decay->Setd0errProngs(3,d0err);
456 // TODO: setter missing
457
458 // update dca for prong combinations
459 Double_t xdummy=0.,ydummy=0.;
460 Double_t dca[3];
461 dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
462 dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
463 dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
464 decay->SetDCAs(3,dca);
465
466 // update dist12 and dist23
467 primaryVertex->GetXYZ(pos);
468 decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
469 +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
470 +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
471 decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
472 +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
473 +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
474
475 delete v123;delete v12;delete v23;
476
477 Double_t px[3],py[3],pz[3];
478 for (Int_t i=0;i<3;++i) {
479 const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
480 px[i]=t->Px();
481 py[i]=t->Py();
482 pz[i]=t->Pz();
483 }
484 decay->SetPxPyPzProngs(3,px,py,pz);
485 }
486 }
487}
488
489void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArray *mcs) {
490 // Early exit, if this track has nothing in common with the ITS
491
492 if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
493 return;
494
495 // Check if the track was already "improved" (this is done with a trick using layer 7 (ie the 8th))
496 if (TESTBIT(track->GetITSClusterMap(),7)) return;
497 //
498
499
500 // Get reconstructed track parameters
501 AliExternalTrackParam et; et.CopyFromVTrack(track);
502 Double_t *param=const_cast<Double_t*>(et.GetParameter());
503//TODO: Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
504
505 // Get MC info
506 Int_t imc=track->GetLabel();
507 if (imc<=0) return;
508 const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
509 Double_t mcx[3];
510 Double_t mcp[3];
511 Double_t mccv[36]={0.};
512 Short_t mcc;
513 mc->XvYvZv(mcx);
514 mc->PxPyPz(mcp);
515 mcc=mc->Charge();
516 AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
517 const Double_t *parammc=mct.GetParameter();
518//TODO: const Double_t *covermc=mct.GetCovariance();
519 AliVertex vtx(mcx,1.,1);
520
521 // Correct reference points and frames according to MC
522 // TODO: B-Field correct?
523 // TODO: failing propagation....
524 et.PropagateToDCA(&vtx,track->GetBz(),10.);
525 et.Rotate(mct.GetAlpha());
526
527 // Select appropriate smearing functions
528 Double_t ptmc=TMath::Abs(mc->Pt());
529 Double_t sd0rpn=0.;
530 Double_t sd0zn =0.;
531 Double_t spt1n =0.;
532 Double_t sd0rpo=0.;
533 Double_t sd0zo =0.;
534 Double_t spt1o =0.;
535 switch (mc->GetPdgCode()) {
536 case 2212: case -2212:
537 sd0rpo=EvalGraph(ptmc,fD0RPResPCur,fD0RPResPCurSA);
538 sd0zo =EvalGraph(ptmc,fD0ZResPCur,fD0ZResPCurSA);
539 spt1o =EvalGraph(ptmc,fPt1ResPCur,fPt1ResPCurSA);
540 sd0rpn=EvalGraph(ptmc,fD0RPResPUpg,fD0RPResPUpgSA);
541 sd0zn =EvalGraph(ptmc,fD0ZResPUpg,fD0ZResPUpgSA);
542 spt1n =EvalGraph(ptmc,fPt1ResPUpg,fPt1ResPUpgSA);
543 break;
544 case 321: case -321:
545 sd0rpo=EvalGraph(ptmc,fD0RPResKCur,fD0RPResKCurSA);
546 sd0zo =EvalGraph(ptmc,fD0ZResKCur,fD0ZResKCurSA);
547 spt1o =EvalGraph(ptmc,fPt1ResKCur,fPt1ResKCurSA);
548 sd0rpn=EvalGraph(ptmc,fD0RPResKUpg,fD0RPResKUpgSA);
549 sd0zn =EvalGraph(ptmc,fD0ZResKUpg,fD0ZResKUpgSA);
550 spt1n =EvalGraph(ptmc,fPt1ResKUpg,fPt1ResKUpgSA);
551 break;
552 case 211: case -211:
553 sd0rpo=EvalGraph(ptmc,fD0RPResPiCur,fD0RPResPiCurSA);
554 sd0zo =EvalGraph(ptmc,fD0ZResPiCur,fD0ZResPiCurSA);
555 spt1o =EvalGraph(ptmc,fPt1ResPiCur,fPt1ResPiCurSA);
556 sd0rpn=EvalGraph(ptmc,fD0RPResPiUpg,fD0RPResPiUpgSA);
557 sd0zn =EvalGraph(ptmc,fD0ZResPiUpg,fD0ZResPiUpgSA);
558 spt1n =EvalGraph(ptmc,fPt1ResPiUpg,fPt1ResPiUpgSA);
559 break;
560 default:
561 return;
562 }
563
564 // Use the same units (i.e. cm and GeV/c)! TODO: pt!
565 sd0rpo*=1.e-4;
566 sd0zo *=1.e-4;
567 sd0rpn*=1.e-4;
568 sd0zn *=1.e-4;
569
570 // Apply the smearing
571 Double_t d0zo =param [1];
572 Double_t d0zmc =parammc[1];
573 Double_t d0rpo =param [0];
574 Double_t d0rpmc=parammc[0];
575 Double_t pt1o =param [4];
576 Double_t pt1mc =parammc[4];
577 Double_t dd0zo =d0zo-d0zmc;
578 Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
579 Double_t d0zn =d0zmc+dd0zn;
580 Double_t dd0rpo=d0rpo-d0rpmc;
581 Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
582 Double_t d0rpn =d0rpmc+dd0rpn;
583 Double_t dpt1o =pt1o-pt1mc;
584 Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
585 Double_t pt1n =pt1mc+dpt1n;
586 param[0]=d0rpn;
587 param[1]=d0zn ;
588 param[4]=pt1n ;
589
590 // Copy the smeared parameters to the AOD track
591 Double_t x[3];
592 Double_t p[3];
593 et.GetXYZ(x);
594 et.GetPxPyPz(p);
595 track->SetPosition(x,kFALSE);
596 track->SetP(p,kTRUE);
597
598
599 // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
600 UChar_t itsClusterMap = track->GetITSClusterMap();
601 SETBIT(itsClusterMap,7);
602 track->SetITSClusterMap(itsClusterMap);
603 //
604
605 // write out debug infos
606 if (fDebugNtuple->GetEntries()<fNDebug) {
607 Int_t idbg=0;
608 fDebugVars[idbg++]=mc->GetPdgCode();
609 fDebugVars[idbg++]=ptmc ;
610 fDebugVars[idbg++]=d0rpo ;
611 fDebugVars[idbg++]=d0zo ;
612 fDebugVars[idbg++]=pt1o ;
613 fDebugVars[idbg++]=sd0rpo;
614 fDebugVars[idbg++]=sd0zo ;
615 fDebugVars[idbg++]=spt1o ;
616 fDebugVars[idbg++]=d0rpn ;
617 fDebugVars[idbg++]=d0zn ;
618 fDebugVars[idbg++]=pt1n ;
619 fDebugVars[idbg++]=sd0rpn;
620 fDebugVars[idbg++]=sd0zn ;
621 fDebugVars[idbg++]=spt1n ;
622 fDebugVars[idbg++]=d0rpmc;
623 fDebugVars[idbg++]=d0zmc ;
624 fDebugVars[idbg++]=pt1mc ;
625 fDebugNtuple->Fill(fDebugVars);
626 PostData(1,fDebugOutput);
627 }
628}
629
630AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
631 //
632 // Helper function to recalculate a vertex.
633 //
634
635 static UShort_t ids[]={1,2,3}; //TODO: unsave...
636 AliVertexerTracks vertexer(bField);
637 vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
638 AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
639 return vertex;
640}
641
642Double_t AliAnalysisTaskSEImproveITS::EvalGraph(Double_t x,const TGraph *graph,const TGraph *graphSA) const {
643 //
644 // Evaluates a TGraph without linear extrapolation. Instead the last
645 // valid point of the graph is used when out of range.
646 // The function assumes an ascending order of X.
647 //
648
649 // TODO: find a pretty solution for this:
650 Int_t n =graph->GetN();
651 Double_t xmin=graph->GetX()[0 ];
652 Double_t xmax=graph->GetX()[n-1];
653 if (x<xmin) {
654 if(!graphSA) return graph->Eval(xmin);
655 Double_t xminSA=graphSA->GetX()[0];
656 if(x<xminSA) return graphSA->Eval(xminSA);
657 return graphSA->Eval(x);
658 }
659 if (x>xmax) return graph->Eval(xmax);
660 return graph->Eval(x);
661}
662
663ClassImp(AliAnalysisTaskSEImproveITS);
664