1 /*************************************************************************
\r
2 * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\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
16 #include <TObjArray.h>
\r
17 #include <TClonesArray.h>
\r
21 #include <TNtuple.h>
\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
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
43 // WARNING: This will affect all tasks in a train after this one
\r
44 // (which is typically desired, though).
\r
47 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS()
\r
48 :AliAnalysisTaskSE(),
\r
85 fRunInVertexing(kFALSE),
\r
86 fImproveTracks(kTRUE),
\r
93 // Default constructor.
\r
97 AliAnalysisTaskSEImproveITS::AliAnalysisTaskSEImproveITS(const char *name,
\r
98 const char *resfileCurURI,
\r
99 const char *resfileUpgURI,
\r
100 Bool_t isRunInVertexing,
\r
102 :AliAnalysisTaskSE(name),
\r
123 fD0ZResPiCurSA (0),
\r
124 fD0RPResPCurSA (0),
\r
125 fD0RPResKCurSA (0),
\r
126 fD0RPResPiCurSA(0),
\r
129 fPt1ResPiCurSA (0),
\r
132 fD0ZResPiUpgSA (0),
\r
133 fD0RPResPUpgSA (0),
\r
134 fD0RPResKUpgSA (0),
\r
135 fD0RPResPiUpgSA(0),
\r
138 fPt1ResPiUpgSA (0),
\r
139 fRunInVertexing(isRunInVertexing),
\r
140 fImproveTracks(kTRUE),
\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
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
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
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
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
231 DefineOutput(1,TNtuple::Class());
\r
234 AliAnalysisTaskSEImproveITS::~AliAnalysisTaskSEImproveITS() {
\r
238 delete fDebugOutput;
\r
241 void AliAnalysisTaskSEImproveITS::UserCreateOutputObjects() {
\r
243 // Creation of user output objects.
\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
250 fDebugOutput->Add(fDebugNtuple );
\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
271 PostData(1,fDebugOutput);
\r
274 void AliAnalysisTaskSEImproveITS::UserExec(Option_t*) {
\r
278 AliAODEvent *ev=0x0;
\r
279 if(!fRunInVertexing) {
\r
280 ev=dynamic_cast<AliAODEvent*>(InputEvent());
\r
282 if(AODEvent() && IsStandardAOD()) ev = dynamic_cast<AliAODEvent*> (AODEvent());
\r
285 Double_t bz=ev->GetMagneticField();
\r
290 // Smear all tracks
\r
291 TClonesArray *mcs=static_cast<TClonesArray*>(ev->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
\r
293 if (fImproveTracks) {
\r
294 for(Int_t itrack=0;itrack<ev->GetNumberOfTracks();++itrack) {
\r
295 SmearTrack(ev->GetTrack(itrack),mcs);
\r
299 // TODO: recalculated primary vertex
\r
300 AliVVertex *primaryVertex=ev->GetPrimaryVertex();
\r
302 // Recalculate all candidates
\r
304 TClonesArray *array2Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("D0toKpi"));
\r
306 for (Int_t icand=0;icand<array2Prong->GetEntries();++icand) {
\r
307 AliAODRecoDecayHF2Prong *decay=static_cast<AliAODRecoDecayHF2Prong*>(array2Prong->At(icand));
\r
309 // recalculate vertices
\r
310 AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
\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
318 ta12.Add(&et1); ta12.Add(&et2);
\r
319 AliESDVertex *v12 =RecalculateVertex(oldSecondaryVertex,&ta12 ,bz);
\r
322 // update secondary vertex
\r
326 decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
\r
327 decay->GetSecondaryVtx()->SetChi2perNDF(v12->GetChi2toNDF());
\r
329 //!!!!TODO: covariance matrix
\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
336 d0err[0] = TMath::Sqrt(covd0z0[0]);
\r
337 et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
\r
339 d0err[1] = TMath::Sqrt(covd0z0[0]);
\r
340 decay->Setd0Prongs(2,d0);
\r
341 decay->Setd0errProngs(2,d0err);
\r
345 Double_t xdummy=0.,ydummy=0.;
\r
347 dca=et1.GetDCA(&et2,bz,xdummy,ydummy);
\r
348 decay->SetDCA(dca);
\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
360 decay->SetPxPyPzProngs(2,px,py,pz);
\r
366 TClonesArray *arrayCascade=static_cast<TClonesArray*>(ev->GetList()->FindObject("Dstar"));
\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
372 AliAODRecoDecayHF2Prong* decay=(AliAODRecoDecayHF2Prong*)decayDstar->Get2Prong();
\r
374 // recalculate vertices
\r
375 //AliVVertex *oldSecondaryVertex=decay->GetSecondaryVtx();
\r
378 AliExternalTrackParam et3; et3.CopyFromVTrack(static_cast<AliAODTrack*>(decayDstar->GetBachelor()));
\r
381 AliNeutralTrackParam *trackD0 = new AliNeutralTrackParam(decay);
\r
383 //!!!!TODO: covariance matrix
\r
386 Double_t d0z0[2],covd0z0[3];
\r
387 Double_t d01[2],d01err[2];
\r
390 et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
\r
392 d01err[0] = TMath::Sqrt(covd0z0[0]);
\r
393 trackD0->PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
\r
395 d01err[1] = TMath::Sqrt(covd0z0[0]);
\r
396 decayDstar->Setd0Prongs(2,d01);
\r
397 decayDstar->Setd0errProngs(2,d01err);
\r
400 delete trackD0; trackD0=NULL;
\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
410 decayDstar->SetPxPyPzProngs(2,px1,py1,pz1);
\r
417 TClonesArray *array3Prong=static_cast<TClonesArray*>(ev->GetList()->FindObject("Charm3Prong"));
\r
419 for (Int_t icand=0;icand<array3Prong->GetEntries();++icand) {
\r
420 AliAODRecoDecayHF3Prong *decay=static_cast<AliAODRecoDecayHF3Prong*>(array3Prong->At(icand));
\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
435 // update secondary vertex
\r
438 decay->GetSecondaryVtx()->SetPosition(pos[0],pos[1],pos[2]);
\r
439 decay->GetSecondaryVtx()->SetChi2perNDF(v123->GetChi2toNDF());
\r
440 //TODO: covariance matrix
\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
447 d0err[0] = TMath::Sqrt(covd0z0[0]);
\r
448 et2.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
\r
450 d0err[1] = TMath::Sqrt(covd0z0[0]);
\r
451 et3.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
\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
458 // update dca for prong combinations
\r
459 Double_t xdummy=0.,ydummy=0.;
\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
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
475 delete v123;delete v12;delete v23;
\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
484 decay->SetPxPyPzProngs(3,px,py,pz);
\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
492 if (!(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)))
\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
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
506 Int_t imc=track->GetLabel();
\r
507 if (imc<=0) return;
\r
508 const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(mcs->At(imc));
\r
511 Double_t mccv[36]={0.};
\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
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
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
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
544 case 321: case -321:
\r
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
552 case 211: case -211:
\r
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
564 // Use the same units (i.e. cm and GeV/c)! TODO: pt!
\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
590 // Copy the smeared parameters to the AOD track
\r
595 track->SetPosition(x,kFALSE);
\r
596 track->SetP(p,kTRUE);
\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
605 // write out debug infos
\r
606 if (fDebugNtuple->GetEntries()<fNDebug) {
\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
630 AliESDVertex* AliAnalysisTaskSEImproveITS::RecalculateVertex(const AliVVertex *old,TObjArray *tracks,Double_t bField) {
\r
632 // Helper function to recalculate a vertex.
\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
642 Double_t AliAnalysisTaskSEImproveITS::EvalGraph(Double_t x,const TGraph *graph,const TGraph *graphSA) const {
\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
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
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
659 if (x>xmax) return graph->Eval(xmax);
\r
660 return graph->Eval(x);
\r
663 ClassImp(AliAnalysisTaskSEImproveITS);
\r