]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSEImproveITS.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added LHC10h run list for flow analysis (Giacomo)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEImproveITS.cxx
CommitLineData
2d4517ec 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>
10371f90 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 "AliAODRecoDecayHF3Prong.h"
2d4517ec 32#include "AliAnalysisTaskSEImproveITS.h"
33
34//
35// Implementation of the "hybrid-approach" for ITS upgrade studies.
36// The tastk smears the track parameters according to estimations
37// from single-track upgrade studies. Afterwards it recalculates
38// the parameters of the reconstructed decays.
39//
40// WARNING: This will affect all tasks in a train after this one
41// (which is typically desired, though).
42//
43
44AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS()
45 :AliAnalysisTaskSE(),
46 fD0ZResPCur (0),
47 fD0ZResKCur (0),
48 fD0ZResPiCur (0),
49 fD0RPResPCur (0),
50 fD0RPResKCur (0),
51 fD0RPResPiCur(0),
52 fPt1ResPCur (0),
53 fPt1ResKCur (0),
54 fPt1ResPiCur (0),
55 fD0ZResPUpg (0),
56 fD0ZResKUpg (0),
57 fD0ZResPiUpg (0),
58 fD0RPResPUpg (0),
59 fD0RPResKUpg (0),
60 fD0RPResPiUpg(0),
61 fPt1ResPUpg (0),
62 fPt1ResKUpg (0),
63 fPt1ResPiUpg (0),
64 fDebugOutput (0),
65 fDebugNtuple (0),
66 fDebugVars (0),
67 fNDebug (0)
68{
69 //
70 // Default constructor.
71 //
72}
73
74AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
75 const char *resfileCurURI,
76 const char *resfileUpgURI,
77 Int_t ndebug)
78 :AliAnalysisTaskSE(name),
79 fD0ZResPCur (0),
80 fD0ZResKCur (0),
81 fD0ZResPiCur (0),
82 fD0RPResPCur (0),
83 fD0RPResKCur (0),
84 fD0RPResPiCur(0),
85 fPt1ResPCur (0),
86 fPt1ResKCur (0),
87 fPt1ResPiCur (0),
88 fD0ZResPUpg (0),
89 fD0ZResKUpg (0),
90 fD0ZResPiUpg (0),
91 fD0RPResPUpg (0),
92 fD0RPResKUpg (0),
93 fD0RPResPiUpg(0),
94 fPt1ResPUpg (0),
95 fPt1ResKUpg (0),
96 fPt1ResPiUpg (0),
97 fDebugOutput (0),
98 fDebugNtuple (0),
99 fDebugVars (0),
100 fNDebug (ndebug)
101{
102 //
103 // Constructor to be used to create the task.
104 // The the URIs specify the resolution files to be used.
105 // They are expected to contain TGraphs with the resolutions
106 // for the current and the upgraded ITS (see code for details).
107 // One may also specify for how many tracks debug information
108 // is written to the output.
109 //
110 TFile *resfileCur=TFile::Open(resfileCurURI);
111 fD0RPResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResP" )));
112 fD0RPResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResK" )));
113 fD0RPResPiCur=new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0RPResPi")));
114 fD0ZResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResP" )));
115 fD0ZResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResK" )));
116 fD0ZResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("D0ZResPi" )));
117 fPt1ResPCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResP" )));
118 fPt1ResKCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResK" )));
119 fPt1ResPiCur =new TGraph(*static_cast<TGraph*>(resfileCur->Get("Pt1ResPi" )));
120 fD0RPResPCur ->SetName("D0RPResPCur" );
121 fD0RPResKCur ->SetName("D0RPResKCur" );
122 fD0RPResPiCur->SetName("D0RPResPiCur");
123 fD0ZResPCur ->SetName("D0ZResPCur" );
124 fD0ZResKCur ->SetName("D0ZResKCur" );
125 fD0ZResPiCur ->SetName("D0ZResPiCur" );
126 fPt1ResPCur ->SetName("Pt1ResPCur" );
127 fPt1ResKCur ->SetName("Pt1ResKCur" );
128 fPt1ResPiCur ->SetName("Pt1ResPiCur" );
129 delete resfileCur;
130 TFile *resfileUpg=TFile::Open(resfileUpgURI );
131 fD0RPResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResP" )));
132 fD0RPResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResK" )));
133 fD0RPResPiUpg=new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0RPResPi")));
134 fD0ZResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResP" )));
135 fD0ZResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResK" )));
136 fD0ZResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("D0ZResPi" )));
137 fPt1ResPUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResP" )));
138 fPt1ResKUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResK" )));
139 fPt1ResPiUpg =new TGraph(*static_cast<TGraph*>(resfileUpg->Get("Pt1ResPi" )));
140 fD0RPResPUpg ->SetName("D0RPResPUpg" );
141 fD0RPResKUpg ->SetName("D0RPResKUpg" );
142 fD0RPResPiUpg->SetName("D0RPResPiUpg");
143 fD0ZResPUpg ->SetName("D0ZResPUpg" );
144 fD0ZResKUpg ->SetName("D0ZResKUpg" );
145 fD0ZResPiUpg ->SetName("D0ZResPiUpg" );
146 fPt1ResPUpg ->SetName("Pt1ResPUpg" );
147 fPt1ResKUpg ->SetName("Pt1ResKUpg" );
148 fPt1ResPiUpg ->SetName("Pt1ResPiUpg" );
149 delete resfileUpg;
150
151 DefineOutput(1,TNtuple::Class());
152}
153
154AliAnalysisTaskSEImproveITS::~AliAnalysisTaskSEImproveITS() {
155 //
156 // Destructor.
157 //
158 delete fDebugOutput;
159}
160
161void AliAnalysisTaskSEImproveITS::UserCreateOutputObjects() {
162 //
163 // Creation of user output objects.
164 //
165 fDebugOutput=new TList();
166 fDebugOutput->SetOwner();
167 fDebugNtuple=new TNtuple("fDebugNtuple","Smearing","pdg:ptmc:d0rpo:d0zo:pt1o:sd0rpo:sd0zo:spt1o:d0rpn:d0zn:pt1n:sd0rpn:sd0zn:spt1n:d0rpmc:d0zmc:pt1mc");
168 fDebugVars=new Float_t[fDebugNtuple->GetNvar()];
169
170 fDebugOutput->Add(fDebugNtuple );
171
172 fDebugOutput->Add(fD0RPResPCur );
173 fDebugOutput->Add(fD0RPResKCur );
174 fDebugOutput->Add(fD0RPResPiCur);
175 fDebugOutput->Add(fD0ZResPCur );
176 fDebugOutput->Add(fD0ZResKCur );
177 fDebugOutput->Add(fD0ZResPiCur );
178 fDebugOutput->Add(fPt1ResPCur );
179 fDebugOutput->Add(fPt1ResKCur );
180 fDebugOutput->Add(fPt1ResPiCur );
181 fDebugOutput->Add(fD0RPResPUpg );
182 fDebugOutput->Add(fD0RPResKUpg );
183 fDebugOutput->Add(fD0RPResPiUpg);
184 fDebugOutput->Add(fD0ZResPUpg );
185 fDebugOutput->Add(fD0ZResKUpg );
186 fDebugOutput->Add(fD0ZResPiUpg );
187 fDebugOutput->Add(fPt1ResPUpg );
188 fDebugOutput->Add(fPt1ResKUpg );
189 fDebugOutput->Add(fPt1ResPiUpg );
190
191 PostData(1,fDebugOutput);
192}
193
194void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
195 //
196 // The event loop
197 //
198 AliAODEvent *ev=dynamic_cast<AliAODEvent*>(InputEvent());
a828d135 199 if(!ev) return;
2d4517ec 200 Double_t bz=ev->GetMagneticField();
201
202 // Smear all tracks
203 TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
204 if (!mcs) return;
205 for (Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack)
206 SmearTrack(ev->GetTrack(itrack),mcs);
207
208 // TODO: recalculated primary vertex
209 AliVVertex *primaryVertex=ev->GetPrimaryVertex();
210
211 // Recalculate all candidates
212 TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
213 if (array3Prong) {
214 for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
215 AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
216
217 // recalculate vertices
218 AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
219 AliExternalTrackParam et1(static_cast<AliAODTrack*>(decay->GetDaughter(0)));
220 AliExternalTrackParam et2(static_cast<AliAODTrack*>(decay->GetDaughter(1)));
221 AliExternalTrackParam et3(static_cast<AliAODTrack*>(decay->GetDaughter(2)));
222 TObjArray ta123,ta12,ta23;
223 ta123.Add(&et1);ta123.Add(&et2);ta123.Add(&et3);
224 ta12. Add(&et1);ta12 .Add(&et2);
225 ta23 .Add(&et2);ta23 .Add(&et3);
226 AliESDVertex *v123=RecalculateVertex(oldSecondaryVertex,&ta123,bz);
227 AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
228 AliESDVertex *v23 =RecalculateVertex(oldSecondaryVertex,&ta23 ,bz);
229
230 // update secondary vertex
231 Double_t pos[3];
232 v123->GetXYZ(pos);
233 decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
234 decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
235 //TODO: covariance matrix
236
237 // update d0 for all progs
238 Double_t d0z0[2],covd0z0[3];
239 Double_t d0[3],d0err[3];
240 et1.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
241 d0[0]=d0z0[0];
242 d0err[0] = TMath::Sqrt(covd0z0[0]);
243 et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
244 d0[1]=d0z0[0];
245 d0err[1] = TMath::Sqrt(covd0z0[0]);
246 et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
247 d0[2]=d0z0[0];
248 d0err[2] = TMath::Sqrt(covd0z0[0]);
249 decay->Setd0Prongs (3,d0 );
250 decay->Setd0errProngs(3,d0err);
251 // TODO: setter missing
252
253 // update dca for prong combinations
254 Double_t xdummy=0.,ydummy=0.;
255 Double_t dca[3];
256 dca[0]=et1.GetDCA(&et2,bz,xdummy,ydummy);
257 dca[1]=et3.GetDCA(&et2,bz,xdummy,ydummy);
258 dca[2]=et1.GetDCA(&et3,bz,xdummy,ydummy);
259 decay->SetDCAs(3,dca);
260
261 // update dist12 and dist23
262 primaryVertex->GetXYZ(pos);
263 decay->SetDist12toPrim(TMath::Sqrt((v12->GetX()-pos[0])*(v12->GetX()-pos[0])
264 +(v12->GetY()-pos[1])*(v12->GetY()-pos[1])
265 +(v12->GetZ()-pos[2])*(v12->GetZ()-pos[2])));
266 decay->SetDist23toPrim(TMath::Sqrt((v23->GetX()-pos[0])*(v23->GetX()-pos[0])
267 +(v23->GetY()-pos[1])*(v23->GetY()-pos[1])
268 +(v23->GetZ()-pos[2])*(v23->GetZ()-pos[2])));
269
270 delete v123;delete v12;delete v23;
271
272 Double_t px[3],py[3],pz[3];
273 for (Int_t i=0;i<3;++i) {
274 const AliAODTrack *t=static_cast<AliAODTrack*>(decay->GetDaughter(i));
275 px[i]=t->Px();
276 py[i]=t->Py();
277 pz[i]=t->Pz();
278 }
279 decay->SetPxPyPzProngs(3,px,py,pz);
280 }
281 }
282}
283
284void AliAnalysisTaskSEImproveITS::SmearTrack(AliAODTrack *track,const TClonesArray *mcs) {
285 // Early exit, if this track has nothing in common with the ITS
286 if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
287 return;
288
289 // Get reconstructed track parameters
290 AliExternalTrackParam et(track);
291 Double_t *param=const_cast<Double_t*>(et.GetParameter());
292//TODO: Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
293
294 // Get MC info
295 Int_t imc=track->GetLabel();
296 if (imc<=0) return;
297 const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
298 Double_t mcx[3];
299 Double_t mcp[3];
300 Double_t mccv[36]={0.};
301 Short_t mcc;
302 mc->XvYvZv(mcx);
303 mc->PxPyPz(mcp);
304 mcc=mc->Charge();
305 AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
306 const Double_t *parammc=mct.GetParameter();
307//TODO: const Double_t *covermc=mct.GetCovariance();
308 AliVertex vtx(mcx,1.,1);
309
310 // Correct reference points and frames according to MC
311 // TODO: B-Field correct?
312 // TODO: failing propagation....
313 et.PropagateToDCA(&vtx,track->GetBz(),10.);
314 et.Rotate(mct.GetAlpha());
315
316 // Select appropriate smearing functions
317 Double_t ptmc=TMath::Abs(mc->Pt());
318 Double_t sd0rpn=0.;
319 Double_t sd0zn =0.;
320 Double_t spt1n =0.;
321 Double_t sd0rpo=0.;
322 Double_t sd0zo =0.;
323 Double_t spt1o =0.;
324 switch (mc->GetPdgCode()) {
325 case 2212: case -2212:
326 sd0rpo=EvalGraph(fD0RPResPCur,ptmc);
327 sd0zo =EvalGraph(fD0ZResPCur ,ptmc);
328 spt1o =EvalGraph(fPt1ResPCur ,ptmc);
329 sd0rpn=EvalGraph(fD0RPResPUpg,ptmc);
330 sd0zn =EvalGraph(fD0ZResPUpg ,ptmc);
331 spt1n =EvalGraph(fPt1ResPUpg ,ptmc);
332 break;
333 case 321: case -321:
334 sd0rpo=EvalGraph(fD0RPResKCur,ptmc);
335 sd0zo =EvalGraph(fD0ZResKCur ,ptmc);
336 spt1o =EvalGraph(fPt1ResKCur ,ptmc);
337 sd0rpn=EvalGraph(fD0RPResKUpg,ptmc);
338 sd0zn =EvalGraph(fD0ZResKUpg ,ptmc);
339 spt1n =EvalGraph(fPt1ResKUpg ,ptmc);
340 break;
341 case 211: case -211:
342 sd0rpo=EvalGraph(fD0RPResPiCur,ptmc);
343 sd0zo =EvalGraph(fD0ZResPiCur ,ptmc);
344 spt1o =EvalGraph(fPt1ResPiCur ,ptmc);
345 sd0rpn=EvalGraph(fD0RPResPiUpg,ptmc);
346 sd0zn =EvalGraph(fD0ZResPiUpg ,ptmc);
347 spt1n =EvalGraph(fPt1ResPiUpg ,ptmc);
348 break;
349 default:
350 return;
351 }
352
353 // Use the same units (i.e. cm and GeV/c)! TODO: pt!
354 sd0rpo*=1.e-4;
355 sd0zo *=1.e-4;
356 sd0rpn*=1.e-4;
357 sd0zn *=1.e-4;
358
359 // Apply the smearing
360 Double_t d0zo =param [1];
361 Double_t d0zmc =parammc[1];
362 Double_t d0rpo =param [0];
363 Double_t d0rpmc=parammc[0];
364 Double_t pt1o =param [4];
365 Double_t pt1mc =parammc[4];
366 Double_t dd0zo =d0zo-d0zmc;
367 Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
368 Double_t d0zn =d0zmc+dd0zn;
369 Double_t dd0rpo=d0rpo-d0rpmc;
370 Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
371 Double_t d0rpn =d0rpmc+dd0rpn;
372 Double_t dpt1o =pt1o-pt1mc;
373 Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
374 Double_t pt1n =pt1mc+dpt1n;
375 param[0]=d0rpn;
376 param[1]=d0zn ;
377 param[4]=pt1n ;
378
379 // Copy the smeared parameters to the AOD track
380 Double_t x[3];
381 Double_t p[3];
382 et.GetXYZ(x);
383 et.GetPxPyPz(p);
384 track->SetPosition(x,kFALSE);
385 track->SetP(p,kTRUE);
386
387 // write out debug infos
388 if (fDebugNtuple->GetEntries()<fNDebug) {
389 Int_t idbg=0;
390 fDebugVars[idbg++]=mc->GetPdgCode();
391 fDebugVars[idbg++]=ptmc ;
392 fDebugVars[idbg++]=d0rpo ;
393 fDebugVars[idbg++]=d0zo ;
394 fDebugVars[idbg++]=pt1o ;
395 fDebugVars[idbg++]=sd0rpo;
396 fDebugVars[idbg++]=sd0zo ;
397 fDebugVars[idbg++]=spt1o ;
398 fDebugVars[idbg++]=d0rpn ;
399 fDebugVars[idbg++]=d0zn ;
400 fDebugVars[idbg++]=pt1n ;
401 fDebugVars[idbg++]=sd0rpn;
402 fDebugVars[idbg++]=sd0zn ;
403 fDebugVars[idbg++]=spt1n ;
404 fDebugVars[idbg++]=d0rpmc;
405 fDebugVars[idbg++]=d0zmc ;
406 fDebugVars[idbg++]=pt1mc ;
407 fDebugNtuple->Fill(fDebugVars);
408 PostData(1,fDebugOutput);
409 }
410}
411
412AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
413 //
414 // Helper function to recalculate a vertex.
415 //
416
417 static UShort_t ids[]={1,2,3}; //TODO: unsave...
418 AliVertexerTracks vertexer(bField);
419 vertexer.SetVtxStart(old->GetX(),old->GetY(),old->GetZ());
420 AliESDVertex *vertex=vertexer.VertexForSelectedTracks(tracks,ids);
421 return vertex;
422}
423
424Double_t AliAnalysisTaskSEImproveITS::EvalGraph(const TGraph *graph,Double_t x) const {
425 //
426 // Evaluates a TGraph without linear extrapolation. Instead the last
427 // valid point of the graph is used when out of range.
428 // The function assumes an ascending order of X.
429 //
430
431 // TODO: find a pretty solution for this:
432 Int_t n =graph->GetN();
433 Double_t xmin=graph->GetX()[0 ];
434 Double_t xmax=graph->GetX()[n-1];
435 if (x<xmin) return graph->Eval(xmin);
436 if (x>xmax) return graph->Eval(xmax);
437 return graph->Eval(x);
438}
439
440ClassImp(AliAnalysisTaskSEImproveITS);
441