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