1 /*************************************************************************
2 * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 #include <TObjArray.h>
17 #include <TClonesArray.h>
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"
32 #include "AliAnalysisTaskSEImproveITS.h"
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.
40 // WARNING: This will affect all tasks in a train after this one
41 // (which is typically desired, though).
44 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS()
70 // Default constructor.
74 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
75 const char *resfileCurURI,
76 const char *resfileUpgURI,
78 :AliAnalysisTaskSE(name),
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.
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" );
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" );
151 DefineOutput(1,TNtuple::Class());
154 AliAnalysisTaskSEImproveITS::~AliAnalysisTaskSEImproveITS() {
161 void AliAnalysisTaskSEImproveITS::UserCreateOutputObjects() {
163 // Creation of user output objects.
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()];
170 fDebugOutput->Add(fDebugNtuple );
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 );
191 PostData(1,fDebugOutput);
194 void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
198 AliAODEvent *ev=dynamic_cast<AliAODEvent*>(InputEvent());
200 Double_t bz=ev->GetMagneticField();
203 TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
205 for (Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack)
206 SmearTrack(ev->GetTrack(itrack),mcs);
208 // TODO: recalculated primary vertex
209 AliVVertex *primaryVertex=ev->GetPrimaryVertex();
211 // Recalculate all candidates
212 TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
214 for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
215 AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
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);
230 // update secondary vertex
233 decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
234 decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
235 //TODO: covariance matrix
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);
242 d0err[0] = TMath::Sqrt(covd0z0[0]);
243 et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
245 d0err[1] = TMath::Sqrt(covd0z0[0]);
246 et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
248 d0err[2] = TMath::Sqrt(covd0z0[0]);
249 decay->Setd0Prongs (3,d0 );
250 decay->Setd0errProngs(3,d0err);
251 // TODO: setter missing
253 // update dca for prong combinations
254 Double_t xdummy=0.,ydummy=0.;
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);
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])));
270 delete v123;delete v12;delete v23;
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));
279 decay->SetPxPyPzProngs(3,px,py,pz);
284 void 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)))
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());
295 Int_t imc=track->GetLabel();
297 const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
300 Double_t mccv[36]={0.};
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);
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());
316 // Select appropriate smearing functions
317 Double_t ptmc=TMath::Abs(mc->Pt());
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);
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);
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);
353 // Use the same units (i.e. cm and GeV/c)! TODO: pt!
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;
379 // Copy the smeared parameters to the AOD track
384 track->SetPosition(x,kFALSE);
385 track->SetP(p,kTRUE);
387 // write out debug infos
388 if (fDebugNtuple->GetEntries()<fNDebug) {
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);
412 AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
414 // Helper function to recalculate a vertex.
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);
424 Double_t AliAnalysisTaskSEImproveITS::EvalGraph(const TGraph *graph,Double_t x) const {
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.
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);
440 ClassImp(AliAnalysisTaskSEImproveITS);