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