]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Coding conventions
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-2008, 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/* $Id$ */
17
18//*************************************************************************
19// Class AliAnalysisTaskSEDplus
20// AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and
21//comparison of heavy-flavour decay candidates
22// to MC truth (kinematics stored in the AOD)
23// Authors: Renu Bala, bala@to.infn.it
24// F. Prino, prino@to.infn.it
25// G. Ortona, ortona@to.infn.it
26/////////////////////////////////////////////////////////////
27
28#include <TClonesArray.h>
29#include <TNtuple.h>
30#include <TCanvas.h>
31#include <TList.h>
32#include <TString.h>
33#include <TDatabasePDG.h>
34
35#include "AliAnalysisManager.h"
36#include "AliRDHFCutsDplustoKpipi.h"
37#include "AliAODHandler.h"
38#include "AliAODEvent.h"
39#include "AliAODVertex.h"
40#include "AliAODTrack.h"
41#include "AliAODRecoDecayHF3Prong.h"
42#include "AliAnalysisVertexingHF.h"
43#include "AliAnalysisTaskSE.h"
44#include "AliAnalysisTaskSEDplus.h"
45#include "AliNormalizationCounter.h"
46
47ClassImp(AliAnalysisTaskSEDplus)
48
49
50//________________________________________________________________________
51AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
52AliAnalysisTaskSE(),
53 fOutput(0),
54 fHistNEvents(0),
55 fPtVsMass(0),
56 fPtVsMassTC(0),
57 fYVsPt(0),
58 fYVsPtTC(0),
59 fYVsPtSig(0),
60 fYVsPtSigTC(0),
61 fNtupleDplus(0),
62 fUpmasslimit(1.965),
63 fLowmasslimit(1.765),
64 fNPtBins(0),
65 fBinWidth(0.002),
66 fListCuts(0),
67 fRDCutsProduction(0),
68 fRDCutsAnalysis(0),
69 fCounter(0),
70 fFillNtuple(kFALSE),
71 fReadMC(kFALSE),
72 fUseStrangeness(kFALSE),
73 fUseBit(kTRUE),
74 fCutsDistr(kFALSE),
75 fDoImpPar(kFALSE),
76 fNImpParBins(400),
77 fLowerImpPar(-2000.),
78 fHigherImpPar(2000.),
79 fDoLS(0)
80{
81 // Default constructor
82}
83
84//________________________________________________________________________
85AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple):
86 AliAnalysisTaskSE(name),
87 fOutput(0),
88 fHistNEvents(0),
89 fPtVsMass(0),
90 fPtVsMassTC(0),
91 fYVsPt(0),
92 fYVsPtTC(0),
93 fYVsPtSig(0),
94 fYVsPtSigTC(0),
95 fNtupleDplus(0),
96 fUpmasslimit(1.965),
97 fLowmasslimit(1.765),
98 fNPtBins(0),
99 fBinWidth(0.002),
100 fListCuts(0),
101 fRDCutsProduction(dpluscutsprod),
102 fRDCutsAnalysis(dpluscutsana),
103 fCounter(0),
104 fFillNtuple(fillNtuple),
105 fReadMC(kFALSE),
106 fUseStrangeness(kFALSE),
107 fUseBit(kTRUE),
108 fCutsDistr(kFALSE),
109 fDoImpPar(kFALSE),
110 fNImpParBins(400),
111 fLowerImpPar(-2000.),
112 fHigherImpPar(2000.),
113 fDoLS(0)
114{
115 //
116 // Standrd constructor
117 //
118 fNPtBins=fRDCutsAnalysis->GetNPtBins();
119 // Default constructor
120 // Output slot #1 writes into a TList container
121 DefineOutput(1,TList::Class()); //My private output
122 // Output slot #2 writes cut to private output
123 // DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
124 DefineOutput(2,TList::Class());
125// Output slot #3 writes cut to private output
126 DefineOutput(3,AliNormalizationCounter::Class());
127
128 if(fFillNtuple){
129 // Output slot #4 writes into a TNtuple container
130 DefineOutput(4,TNtuple::Class()); //My private output
131 }
132}
133
134//________________________________________________________________________
135AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
136{
137 //
138 // Destructor
139 //
140 if (fOutput) {
141 delete fOutput;
142 fOutput = 0;
143 }
144 if(fHistNEvents){
145 delete fHistNEvents;
146 fHistNEvents=0;
147 }
148 for(Int_t i=0;i<3;i++){
149 if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;}
150 }
151
152 for(Int_t i=0;i<3*fNPtBins;i++){
153 if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
154 if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
155 if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
156 if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
157 if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
158 if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
159 if(fPtKHist[i]){ delete fPtKHist[i]; fPtKHist[i]=0;}
160 if(fPtpi1Hist[i]){ delete fPtpi1Hist[i]; fPtpi1Hist[i]=0;}
161 if(fPtpi2Hist[i]){ delete fPtpi2Hist[i]; fPtpi2Hist[i]=0;}
162 if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
163 if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
164 if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
165 if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
166
167 if(fDLxy[i]){delete fDLxy[i]; fDLxy[i]=0;}
168 if(fDLxyTC[i]){delete fDLxyTC[i]; fDLxyTC[i]=0;}
169 if(fCosxy[i]){delete fCosxy[i]; fCosxy[i]=0;}
170 if(fCosxyTC[i]){delete fCosxyTC[i]; fCosxyTC[i]=0;}
171 if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
172 if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
173 if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
174 if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
175 if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
176 if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
177 if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
178 if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
179 }
180
181 for(Int_t i=0;i<3;i++){
182 if(fCorreld0Kd0pi[i]){ delete fCorreld0Kd0pi[i]; fCorreld0Kd0pi[i]=0;}
183 }
184
185 if(fPtVsMass){
186 delete fPtVsMass;
187 fPtVsMass=0;
188 }
189 if(fPtVsMassTC){
190 delete fPtVsMassTC;
191 fPtVsMassTC=0;
192 }
193 if(fYVsPt){
194 delete fYVsPt;
195 fYVsPt=0;
196 }
197 if(fYVsPtTC){
198 delete fYVsPtTC;
199 fYVsPtTC=0;
200 }
201 if(fYVsPtSig){
202 delete fYVsPtSig;
203 fYVsPtSig=0;
204 }
205 if(fYVsPtSigTC){
206 delete fYVsPtSigTC;
207 fYVsPtSigTC=0;
208 }
209 if(fNtupleDplus){
210 delete fNtupleDplus;
211 fNtupleDplus=0;
212 }
213 if (fListCuts) {
214 delete fListCuts;
215 fListCuts = 0;
216 }
217 if(fRDCutsProduction){
218 delete fRDCutsProduction;
219 fRDCutsProduction = 0;
220 }
221 if(fRDCutsAnalysis){
222 delete fRDCutsAnalysis;
223 fRDCutsAnalysis = 0;
224 }
225 for(Int_t i=0; i<5; i++){
226 delete fHistMassPtImpParTC[i];
227 }
228 if(fCounter){
229 delete fCounter;
230 fCounter = 0;
231 }
232
233
234}
235//_________________________________________________________________
236void AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
237 // set invariant mass limits
238 Float_t bw=GetBinWidth();
239 fUpmasslimit = 1.865+range;
240 fLowmasslimit = 1.865-range;
241 SetBinWidth(bw);
242}
243//_________________________________________________________________
244void AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
245 // set invariant mass limits
246 if(uplimit>lowlimit)
247 {
248 Float_t bw=GetBinWidth();
249 fUpmasslimit = lowlimit;
250 fLowmasslimit = uplimit;
251 SetBinWidth(bw);
252 }
253}
254//________________________________________________________________
255void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
256 Float_t width=w;
257 Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
258 Int_t missingbins=4-nbins%4;
259 nbins=nbins+missingbins;
260 width=(fUpmasslimit-fLowmasslimit)/nbins;
261 if(missingbins!=0){
262 printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
263 }
264 else{
265 if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
266 }
267 fBinWidth=width;
268}
269//_________________________________________________________________
270Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
271 return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
272}
273//_________________________________________________________________
274void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
275 //
276 //
277 // Fill the Like Sign histograms
278 //
279 if(fDebug>1)printf("started LS\n");
280
281 //histograms for like sign
282 Int_t nbins=GetNBinsHistos();;
283 TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
284 TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
285
286 Int_t nPosTrks=0,nNegTrks=0;
287 Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
288 Int_t nDplusLS=0;
289 Int_t nDminusLS=0;
290 Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
291 Int_t index=0;
292
293 // loop over like sign candidates
294 for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
295 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
296 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
297 Bool_t unsetvtx=kFALSE;
298 if(!d->GetOwnPrimaryVtx()) {
299 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
300 unsetvtx=kTRUE;
301 }
302 Double_t ptCand = d->Pt();
303 Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
304 if(iPtBin<0)continue;
305 Int_t sign= d->GetCharge();
306 if(sign>0){
307 nPosTrks++;
308 }else{
309 nNegTrks++;
310 }
311 // if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){
312 fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
313 Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
314 if(passTightCuts>0){
315 Float_t invMass = d->InvMassDplus();
316 index=GetLSHistoIndex(iPtBin);
317 fMassHistLS[index+1]->Fill(invMass);
318 if(sign>0){
319 histLSPlus->Fill(invMass);
320 nDplusLS++;
321 }else{
322 histLSMinus->Fill(invMass);
323 nDminusLS++;
324 }
325 if(fCutsDistr){
326 Double_t dlen=d->DecayLength();
327 Double_t cosp=d->CosPointingAngle();
328 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
329 Double_t dca=d->GetDCA();
330 Double_t sigvert=d->GetSigmaVert();
331 Double_t ptmax=0;
332 for(Int_t i=0;i<3;i++){
333 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
334 }
335 fCosPHistLS[iPtBin]->Fill(cosp);
336 fDLenHistLS[iPtBin]->Fill(dlen);
337 fSumd02HistLS[iPtBin]->Fill(sumD02);
338 fSigVertHistLS[iPtBin]->Fill(sigvert);
339 fPtMaxHistLS[iPtBin]->Fill(ptmax);
340 fDCAHistLS[iPtBin]->Fill(dca);
341 }
342 }
343 if(unsetvtx) d->UnsetOwnPrimaryVtx();
344 }
345 //wei:normalized to the number of combinations (production)
346 //wei2:normalized to the number of LS/OS (production)
347 //wei3:normalized to the number of LS/OS (analysis)
348 //wei4:normalized to the number of combinations (analysis)
349 Float_t wei2=0;
350 if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
351 Float_t wei3=0;
352 if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)(nDplusLS+nDminusLS);
353 Float_t weiplus=1.,weiminus=1.;
354 Float_t wei4plus=1.,wei4minus=1.;
355 //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
356 if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
357 if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
358 if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
359 if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);
360
361 fMassHistLS[index]->Add(histLSPlus,weiplus);
362 fMassHistLS[index]->Add(histLSMinus,weiminus);
363 fMassHistLS[index+2]->Add(histLSPlus,wei2);
364 fMassHistLS[index+2]->Add(histLSMinus,wei2);
365 fMassHistLS[index+3]->Add(histLSPlus,wei3);
366 fMassHistLS[index+3]->Add(histLSMinus,wei3);
367 fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
368 fMassHistLS[index+4]->Add(histLSMinus,wei4minus);
369
370 delete histLSPlus;histLSPlus=0;
371 delete histLSMinus;histLSMinus=0;
372
373 if(fDebug>1) printf("LS analysis terminated\n");
374}
375
376//__________________________________________
377void AliAnalysisTaskSEDplus::Init(){
378 //
379 // Initialization
380 //
381 if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
382
383 //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
384 fListCuts=new TList();
385 AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction);
386 production->SetName("ProductionCuts");
387 AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
388 analysis->SetName("AnalysisCuts");
389
390 fListCuts->Add(production);
391 fListCuts->Add(analysis);
392 PostData(2,fListCuts);
393
394 return;
395}
396
397//________________________________________________________________________
398void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
399{
400 // Create the output container
401 //
402 if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
403
404 // Several histograms are more conveniently managed in a TList
405 fOutput = new TList();
406 fOutput->SetOwner();
407 fOutput->SetName("OutputHistos");
408
409 TString hisname;
410 Int_t index=0;
411 Int_t nbins=GetNBinsHistos();
412 fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5);
413 fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5);
414 fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5);
415 for(Int_t i=0;i<3;i++){
416 fHistCentrality[i]->Sumw2();
417 fOutput->Add(fHistCentrality[i]);
418 }
419 for(Int_t i=0;i<fNPtBins;i++){
420
421 index=GetHistoIndex(i);
422
423 hisname.Form("hMassPt%d",i);
424 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
425 fMassHist[index]->Sumw2();
426 hisname.Form("hCosPAllPt%d",i);
427 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
428 fCosPHist[index]->Sumw2();
429 hisname.Form("hDLenAllPt%d",i);
430 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
431 fDLenHist[index]->Sumw2();
432 hisname.Form("hSumd02AllPt%d",i);
433 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
434 fSumd02Hist[index]->Sumw2();
435 hisname.Form("hSigVertAllPt%d",i);
436 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
437 fSigVertHist[index]->Sumw2();
438 hisname.Form("hPtMaxAllPt%d",i);
439 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
440 fPtMaxHist[index]->Sumw2();
441 hisname.Form("hPtKPt%d",i);
442 fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
443 fPtKHist[index]->Sumw2();
444 hisname.Form("hPtpi1Pt%d",i);
445 fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
446 fPtpi1Hist[index]->Sumw2();
447 hisname.Form("hPtpi2Pt%d",i);
448 fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
449 fPtpi2Hist[index]->Sumw2();
450 hisname.Form("hDCAAllPt%d",i);
451 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
452 fDCAHist[index]->Sumw2();
453
454 hisname.Form("hDLxyPt%d",i);
455 fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
456 fDLxy[index]->Sumw2();
457 hisname.Form("hCosxyPt%d",i);
458 fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
459 fCosxy[index]->Sumw2();
460 hisname.Form("hDLxyPt%dTC",i);
461 fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
462 fDLxyTC[index]->Sumw2();
463 hisname.Form("hCosxyPt%dTC",i);
464 fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
465 fCosxyTC[index]->Sumw2();
466
467 hisname.Form("hMassPt%dTC",i);
468 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
469 fMassHistTC[index]->Sumw2();
470 hisname.Form("hMassPt%dTCPlus",i);
471 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
472 fMassHistTCPlus[index]->Sumw2();
473 hisname.Form("hMassPt%dTCMinus",i);
474 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
475 fMassHistTCMinus[index]->Sumw2();
476
477
478
479 index=GetSignalHistoIndex(i);
480 hisname.Form("hSigPt%d",i);
481 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
482 fMassHist[index]->Sumw2();
483 hisname.Form("hCosPSigPt%d",i);
484 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
485 fCosPHist[index]->Sumw2();
486 hisname.Form("hDLenSigPt%d",i);
487 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
488 fDLenHist[index]->Sumw2();
489 hisname.Form("hSumd02SigPt%d",i);
490 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
491 fSumd02Hist[index]->Sumw2();
492 hisname.Form("hSigVertSigPt%d",i);
493 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
494 fSigVertHist[index]->Sumw2();
495 hisname.Form("hPtMaxSigPt%d",i);
496 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
497 fPtMaxHist[index]->Sumw2();
498 hisname.Form("hPtKSigPt%d",i);
499 fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
500 fPtKHist[index]->Sumw2();
501 hisname.Form("hPtpi1SigPt%d",i);
502 fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
503 fPtpi1Hist[index]->Sumw2();
504 hisname.Form("hPtpi2SigPt%d",i);
505 fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
506 fPtpi2Hist[index]->Sumw2();
507
508 hisname.Form("hDCASigPt%d",i);
509 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
510 fDCAHist[index]->Sumw2();
511
512 hisname.Form("hDLxySigPt%d",i);
513 fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
514 fDLxy[index]->Sumw2();
515 hisname.Form("hCosxySigPt%d",i);
516 fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
517 fCosxy[index]->Sumw2();
518 hisname.Form("hDLxySigPt%dTC",i);
519 fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
520 fDLxyTC[index]->Sumw2();
521 hisname.Form("hCosxySigPt%dTC",i);
522 fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
523 fCosxyTC[index]->Sumw2();
524 hisname.Form("hSigPt%dTC",i);
525 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
526 fMassHistTC[index]->Sumw2();
527 hisname.Form("hSigPt%dTCPlus",i);
528 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
529 fMassHistTCPlus[index]->Sumw2();
530 hisname.Form("hSigPt%dTCMinus",i);
531 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
532 fMassHistTCMinus[index]->Sumw2();
533
534
535 index=GetBackgroundHistoIndex(i);
536 hisname.Form("hBkgPt%d",i);
537 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
538 fMassHist[index]->Sumw2();
539 hisname.Form("hCosPBkgPt%d",i);
540 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
541 fCosPHist[index]->Sumw2();
542 hisname.Form("hDLenBkgPt%d",i);
543 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
544 fDLenHist[index]->Sumw2();
545 hisname.Form("hSumd02BkgPt%d",i);
546 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
547 fSumd02Hist[index]->Sumw2();
548 hisname.Form("hSigVertBkgPt%d",i);
549 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
550 fSigVertHist[index]->Sumw2();
551 hisname.Form("hPtMaxBkgPt%d",i);
552 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
553 fPtMaxHist[index]->Sumw2();
554 hisname.Form("hPtKBkgPt%d",i);
555 fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
556 fPtKHist[index]->Sumw2();
557 hisname.Form("hPtpi1BkgPt%d",i);
558 fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
559 fPtpi1Hist[index]->Sumw2();
560 hisname.Form("hPtpi2BkgPt%d",i);
561 fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
562 fPtpi2Hist[index]->Sumw2();
563 hisname.Form("hDCABkgPt%d",i);
564 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
565 fDCAHist[index]->Sumw2();
566
567 hisname.Form("hDLxyBkgPt%d",i);
568 fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
569 fDLxy[index]->Sumw2();
570 hisname.Form("hCosxyBkgPt%d",i);
571 fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
572 fCosxy[index]->Sumw2();
573 hisname.Form("hDLxyBkgPt%dTC",i);
574 fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
575 fDLxyTC[index]->Sumw2();
576 hisname.Form("hCosxyBkgPt%dTC",i);
577 fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
578 fCosxyTC[index]->Sumw2();
579
580
581 hisname.Form("hBkgPt%dTC",i);
582 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
583 fMassHistTC[index]->Sumw2();
584 hisname.Form("hBkgPt%dTCPlus",i);
585 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
586 fMassHistTCPlus[index]->Sumw2();
587 hisname.Form("hBkgPt%dTCMinus",i);
588 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
589 fMassHistTCMinus[index]->Sumw2();
590 }
591
592
593 for(Int_t i=0; i<3*fNPtBins; i++){
594 fOutput->Add(fMassHist[i]);
595 if(fCutsDistr){
596 fOutput->Add(fCosPHist[i]);
597 fOutput->Add(fDLenHist[i]);
598 fOutput->Add(fSumd02Hist[i]);
599 fOutput->Add(fSigVertHist[i]);
600 fOutput->Add(fPtMaxHist[i]);
601 fOutput->Add(fPtKHist[i]);
602 fOutput->Add(fPtpi1Hist[i]);
603 fOutput->Add(fPtpi2Hist[i]);
604 fOutput->Add(fDCAHist[i]);
605 fOutput->Add(fDLxy[i]);
606 fOutput->Add(fDLxyTC[i]);
607 fOutput->Add(fCosxy[i]);
608 fOutput->Add(fCosxyTC[i]);
609 }
610 fOutput->Add(fMassHistTC[i]);
611 fOutput->Add(fMassHistTCPlus[i]);
612 fOutput->Add(fMassHistTCMinus[i]);
613
614 }
615
616 if(fCutsDistr){
617 fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
618 fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
619 fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
620 for(Int_t i=0; i<3; i++){
621 fCorreld0Kd0pi[i]->Sumw2();
622 fOutput->Add(fCorreld0Kd0pi[i]);
623 }
624 }
625
626 fHistNEvents = new TH1F("fHistNEvents", "number of events ",9,-0.5,8.5);
627 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
628 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
629 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger");
630 fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events");
631 fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate");
632 fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts");
633 fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts");
634 fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
635 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
636 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
637 fHistNEvents->Sumw2();
638 fHistNEvents->SetMinimum(0);
639 fOutput->Add(fHistNEvents);
640
641 fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
642 fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
643 fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
644 fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
645 fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
646 fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
647
648 fOutput->Add(fPtVsMass);
649 fOutput->Add(fPtVsMassTC);
650 fOutput->Add(fYVsPt);
651 fOutput->Add(fYVsPtTC);
652 fOutput->Add(fYVsPtSig);
653 fOutput->Add(fYVsPtSigTC);
654
655
656 //Counter for Normalization
657 TString normName="NormalizationCounter";
658 AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
659 if(cont)normName=(TString)cont->GetName();
660 fCounter = new AliNormalizationCounter(normName.Data());
661 fCounter->Init();
662
663 if(fDoLS) CreateLikeSignHistos();
664 if(fDoImpPar) CreateImpactParameterHistos();
665
666 if(fFillNtuple){
667 OpenFile(4); // 4 is the slot number of the ntuple
668
669 fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
670
671 }
672
673 return;
674}
675
676//________________________________________________________________________
677void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
678{
679 // Execute analysis for current event:
680 // heavy flavor candidates association to MC truth
681
682 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
683
684
685
686 TClonesArray *array3Prong = 0;
687 TClonesArray *arrayLikeSign =0;
688 if(!aod && AODEvent() && IsStandardAOD()) {
689 // In case there is an AOD handler writing a standard AOD, use the AOD
690 // event in memory rather than the input (ESD) event.
691 aod = dynamic_cast<AliAODEvent*> (AODEvent());
692 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
693 // have to taken from the AOD event hold by the AliAODExtension
694 AliAODHandler* aodHandler = (AliAODHandler*)
695 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
696 if(aodHandler->GetExtensions()) {
697 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
698 AliAODEvent *aodFromExt = ext->GetAOD();
699 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
700 arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
701 }
702 } else if(aod) {
703 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
704 arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
705 }
706
707 if(!aod || !array3Prong) {
708 printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
709 return;
710 }
711 if(!arrayLikeSign && fDoLS) {
712 printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
713 return;
714 }
715
716
717 // fix for temporary bug in ESDfilter
718 // the AODs with null vertex pointer didn't pass the PhysSel
719 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
720 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
721 fHistNEvents->Fill(0); // count event
722
723 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
724 Bool_t isEvSelP=kTRUE;
725 isEvSelP=fRDCutsProduction->IsEventSelected(aod); // to have proper PID object settings
726
727 Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
728 fHistCentrality[0]->Fill(centrality);
729 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
730 TString trigclass=aod->GetFiredTriggerClasses();
731 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
732 if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3);
733 if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
734
735 // Post the data already here
736 PostData(1,fOutput);
737 if(!isEvSel)return;
738
739 fHistCentrality[1]->Fill(centrality);
740 fHistNEvents->Fill(1);
741
742 TClonesArray *arrayMC=0;
743 AliAODMCHeader *mcHeader=0;
744
745 // AOD primary vertex
746 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
747 // vtx1->Print();
748 TString primTitle = vtx1->GetTitle();
749 //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
750
751 // load MC particles
752 if(fReadMC){
753
754 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
755 if(!arrayMC) {
756 printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
757 return;
758 }
759
760 // load MC header
761 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
762 if(!mcHeader) {
763 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
764 return;
765 }
766 }
767
768 Int_t n3Prong = array3Prong->GetEntriesFast();
769 // printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
770
771
772 Int_t nOS=0;
773 Int_t index;
774 Int_t pdgDgDplustoKpipi[3]={321,211,211};
775
776 if(fDoLS>1){
777 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
778 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
779 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
780 if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod))nOS++;
781 }
782 }
783 }else{
784 // Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};//TO REMOVE
785 //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
786 Int_t nSelectedloose=0,nSelectedtight=0;
787 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
788 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
789 fHistNEvents->Fill(4);
790 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
791 fHistNEvents->Fill(8);
792 continue;
793 }
794 Bool_t unsetvtx=kFALSE;
795 if(!d->GetOwnPrimaryVtx()){
796 d->SetOwnPrimaryVtx(vtx1);
797 unsetvtx=kTRUE;
798 }
799
800 if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
801
802 Double_t ptCand = d->Pt();
803 Int_t iPtBin = fRDCutsProduction->PtBin(ptCand);
804
805 Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
806 Bool_t recVtx=kFALSE;
807 AliAODVertex *origownvtx=0x0;
808 if(fRDCutsProduction->GetIsPrimaryWithoutDaughters()){
809 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
810 if(fRDCutsProduction->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
811 else fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
812 }
813
814 Int_t labDp=-1;
815 Bool_t isPrimary=kTRUE;
816 Float_t deltaPx=0.;
817 Float_t deltaPy=0.;
818 Float_t deltaPz=0.;
819 Float_t truePt=0.;
820 Float_t xDecay=0.;
821 Float_t yDecay=0.;
822 Float_t zDecay=0.;
823 Float_t pdgCode=-2;
824 Float_t trueImpParXY=0.;
825 if(fReadMC){
826 labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
827 if(labDp>=0){
828 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
829 if(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE;
830 AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
831 deltaPx=partDp->Px()-d->Px();
832 deltaPy=partDp->Py()-d->Py();
833 deltaPz=partDp->Pz()-d->Pz();
834 truePt=partDp->Pt();
835 xDecay=dg0->Xv();
836 yDecay=dg0->Yv();
837 zDecay=dg0->Zv();
838 pdgCode=TMath::Abs(partDp->GetPdgCode());
839 if(!isPrimary){
840 trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
841 }
842 }else{
843 pdgCode=-1;
844 }
845 }
846
847 Double_t invMass=d->InvMassDplus();
848 Double_t rapid=d->YDplus();
849 fYVsPt->Fill(ptCand,rapid);
850 if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid);nOS++;}
851 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
852 if(isFidAcc){
853 fPtVsMass->Fill(invMass,ptCand);
854 if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
855 }
856 Float_t tmp[24];
857 Double_t dlen=d->DecayLength();
858 Double_t cosp=d->CosPointingAngle();
859 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
860 Double_t dca=d->GetDCA();
861 Double_t sigvert=d->GetSigmaVert();
862 Double_t ptmax=0;
863 for(Int_t i=0;i<3;i++){
864 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
865 }
866 Double_t impparXY=d->ImpParXY()*10000.;
867 Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
868 Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
869 if(fFillNtuple){
870 tmp[0]=pdgCode;
871 tmp[1]=deltaPx;
872 tmp[2]=deltaPy;
873 tmp[3]=deltaPz;
874 tmp[4]=truePt;
875 tmp[5]=xDecay;
876 tmp[6]=yDecay;
877 tmp[7]=zDecay;
878 tmp[8]=d->PtProng(0);
879 tmp[9]=d->PtProng(1);
880 tmp[10]=d->PtProng(2);
881 tmp[11]=d->Pt();
882 tmp[12]=cosp;
883 tmp[13]=dlen;
884 tmp[14]=d->Xv();
885 tmp[15]=d->Yv();
886 tmp[16]=d->Zv();
887 tmp[17]=d->InvMassDplus();
888 tmp[18]=sigvert;
889 tmp[19]=d->Getd0Prong(0);
890 tmp[20]=d->Getd0Prong(1);
891 tmp[21]=d->Getd0Prong(2);
892 tmp[22]=dca;
893 tmp[23]=d->Prodd0d0();
894 fNtupleDplus->Fill(tmp);
895 PostData(4,fNtupleDplus);
896 }
897 if(iPtBin>=0){
898 Float_t dlxy=d->NormalizedDecayLengthXY();
899 Float_t cxy=d->CosPointingAngleXY();
900 index=GetHistoIndex(iPtBin);
901 if(isFidAcc){
902 fHistNEvents->Fill(5);
903 nSelectedloose++;
904 fMassHist[index]->Fill(invMass);
905 if(fCutsDistr){
906 fCosPHist[index]->Fill(cosp);
907 fDLenHist[index]->Fill(dlen);
908 fSumd02Hist[index]->Fill(sumD02);
909 fSigVertHist[index]->Fill(sigvert);
910 fPtMaxHist[index]->Fill(ptmax);
911 fPtKHist[index]->Fill(d->PtProng(1));
912 fPtpi1Hist[index]->Fill(d->PtProng(0));
913 fPtpi2Hist[index]->Fill(d->PtProng(2));
914 fDCAHist[index]->Fill(dca);
915 fDLxy[index]->Fill(dlxy);
916 fCosxy[index]->Fill(cxy);
917 fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
918 d->Getd0Prong(2)*d->Getd0Prong(1));
919 }
920 if(passTightCuts){ fHistNEvents->Fill(6);
921 nSelectedtight++;
922 fMassHistTC[index]->Fill(invMass);
923 if(fCutsDistr){
924 fDLxyTC[index]->Fill(dlxy);
925 fCosxyTC[index]->Fill(cxy);
926 }
927 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
928 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
929 if(fDoImpPar){
930 fHistMassPtImpParTC[0]->Fill(arrayForSparse);
931 }
932 }
933 }
934
935 if(fReadMC){
936 // if(fCutsDistr){
937 if(labDp>=0) {
938 index=GetSignalHistoIndex(iPtBin);
939 if(isFidAcc){
940 Float_t factor[3]={1.,1.,1.};
941 if(fUseStrangeness){
942 for(Int_t iprong=0;iprong<3;iprong++){
943 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
944 Int_t labd= trad->GetLabel();
945 if(labd>=0){
946 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
947 if(dau){
948 Int_t labm = dau->GetMother();
949 if(labm>=0){
950 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
951 if(mot){
952 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
953 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
954 else factor[iprong]=1./.6;
955 // fNentries->Fill(12);
956 }
957 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
958 factor[iprong]=1./0.25;
959 // fNentries->Fill(13);
960 }//if 3122
961 }//if(mot)
962 }//if labm>0
963 }//if(dau)
964 }//if labd>=0
965 }//prong loop
966 }
967 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
968 fMassHist[index]->Fill(invMass);
969 if(fCutsDistr){
970 fCosPHist[index]->Fill(cosp,fact);
971 fDLenHist[index]->Fill(dlen,fact);
972 fDLxy[index]->Fill(dlxy);
973 fCosxy[index]->Fill(cxy);
974
975 Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2];
976 fSumd02Hist[index]->Fill(sumd02s);
977 fSigVertHist[index]->Fill(sigvert,fact);
978 fPtMaxHist[index]->Fill(ptmax,fact);
979 fPtKHist[index]->Fill(d->PtProng(1),fact);
980 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
981 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
982 fDCAHist[index]->Fill(dca,fact);
983 fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
984 d->Getd0Prong(2)*d->Getd0Prong(1));
985 }
986 if(passTightCuts){
987 fMassHistTC[index]->Fill(invMass);
988 if(fCutsDistr){
989 fDLxyTC[index]->Fill(dlxy);
990 fCosxyTC[index]->Fill(cxy);
991 }
992 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
993 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
994 if(fDoImpPar){
995 if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
996 else{
997 fHistMassPtImpParTC[2]->Fill(arrayForSparse);
998 fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
999 }
1000 }
1001 }
1002 }
1003 fYVsPtSig->Fill(ptCand,rapid);
1004 if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
1005 }else{
1006 index=GetBackgroundHistoIndex(iPtBin);
1007 if(isFidAcc){
1008 Float_t factor[3]={1.,1.,1.};
1009 if(fUseStrangeness){
1010 for(Int_t iprong=0;iprong<3;iprong++){
1011 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1012 Int_t labd= trad->GetLabel();
1013 if(labd>=0){
1014 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1015 if(dau){
1016 Int_t labm = dau->GetMother();
1017 if(labm>=0){
1018 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1019 if(mot){
1020 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1021 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1022 else factor[iprong]=1./.6;
1023 // fNentries->Fill(12);
1024 }
1025 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1026 factor[iprong]=1./0.25;
1027 // fNentries->Fill(13);
1028 }//if 3122
1029 }//if(mot)
1030 }//if labm>0
1031 }//if(dau)
1032 }//if labd>=0
1033 }//prong loop
1034 }
1035
1036 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1037 fMassHist[index]->Fill(invMass);
1038 if(fCutsDistr){
1039 fCosPHist[index]->Fill(cosp,fact);
1040 fDLenHist[index]->Fill(dlen,fact);
1041 fDLxy[index]->Fill(dlxy);
1042 fCosxy[index]->Fill(cxy);
1043
1044 Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2];
1045 fSumd02Hist[index]->Fill(sumd02s);
1046 fSigVertHist[index]->Fill(sigvert,fact);
1047 fPtMaxHist[index]->Fill(ptmax,fact);
1048 fPtKHist[index]->Fill(d->PtProng(1),fact);
1049 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1050 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1051 fDCAHist[index]->Fill(dca,fact);
1052 fCorreld0Kd0pi[2]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1053 d->Getd0Prong(2)*d->Getd0Prong(1));
1054 }
1055 if(passTightCuts){
1056 fMassHistTC[index]->Fill(invMass);
1057 if(fCutsDistr){
1058 fDLxyTC[index]->Fill(dlxy);
1059 fCosxyTC[index]->Fill(cxy);
1060 }
1061 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1062 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1063 if(fDoImpPar){
1064 fHistMassPtImpParTC[4]->Fill(arrayForSparse);
1065 }
1066 }
1067 }
1068 }
1069
1070 }
1071 }
1072
1073 if(recVtx)fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
1074 }
1075 if(unsetvtx) d->UnsetOwnPrimaryVtx();
1076 }
1077 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1078 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1079 }
1080 //start LS analysis
1081 if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1082
1083 PostData(1,fOutput);
1084 PostData(2,fListCuts);
1085 PostData(3,fCounter);
1086 return;
1087}
1088
1089//________________________________________________________________________
1090void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){
1091 // Histos for Like Sign bckground
1092
1093 TString hisname;
1094 Int_t indexLS=0;
1095 Int_t index=0;
1096 Int_t nbins=GetNBinsHistos();
1097 for(Int_t i=0;i<fNPtBins;i++){
1098
1099 index=GetHistoIndex(i);
1100 indexLS=GetLSHistoIndex(i);
1101
1102 hisname.Form("hLSPt%dLC",i);
1103 fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1104 fMassHistLS[indexLS]->Sumw2();
1105 hisname.Form("hLSPt%dTC",i);
1106 fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1107 fMassHistLSTC[indexLS]->Sumw2();
1108
1109 hisname.Form("hCosPAllPt%dLS",i);
1110 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1111 fCosPHistLS[index]->Sumw2();
1112 hisname.Form("hDLenAllPt%dLS",i);
1113 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1114 fDLenHistLS[index]->Sumw2();
1115 hisname.Form("hSumd02AllPt%dLS",i);
1116 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1117 fSumd02HistLS[index]->Sumw2();
1118 hisname.Form("hSigVertAllPt%dLS",i);
1119 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1120 fSigVertHistLS[index]->Sumw2();
1121 hisname.Form("hPtMaxAllPt%dLS",i);
1122 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1123 fPtMaxHistLS[index]->Sumw2();
1124 hisname.Form("hDCAAllPt%dLS",i);
1125 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1126 fDCAHistLS[index]->Sumw2();
1127
1128 index=GetSignalHistoIndex(i);
1129 indexLS++;
1130
1131 hisname.Form("hLSPt%dLCnw",i);
1132 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1133 fMassHistLS[indexLS]->Sumw2();
1134 hisname.Form("hLSPt%dTCnw",i);
1135 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1136 fMassHistLSTC[indexLS]->Sumw2();
1137
1138 hisname.Form("hCosPSigPt%dLS",i);
1139 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1140 fCosPHistLS[index]->Sumw2();
1141 hisname.Form("hDLenSigPt%dLS",i);
1142 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1143 fDLenHistLS[index]->Sumw2();
1144 hisname.Form("hSumd02SigPt%dLS",i);
1145 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1146 fSumd02HistLS[index]->Sumw2();
1147 hisname.Form("hSigVertSigPt%dLS",i);
1148 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1149 fSigVertHistLS[index]->Sumw2();
1150 hisname.Form("hPtMaxSigPt%dLS",i);
1151 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1152 fPtMaxHistLS[index]->Sumw2();
1153 hisname.Form("hDCASigPt%dLS",i);
1154 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1155 fDCAHistLS[index]->Sumw2();
1156
1157 index=GetBackgroundHistoIndex(i);
1158 indexLS++;
1159
1160 hisname.Form("hLSPt%dLCntrip",i);
1161 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1162 fMassHistLS[indexLS]->Sumw2();
1163 hisname.Form("hLSPt%dTCntrip",i);
1164 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1165 fMassHistLSTC[indexLS]->Sumw2();
1166
1167 hisname.Form("hCosPBkgPt%dLS",i);
1168 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1169 fCosPHistLS[index]->Sumw2();
1170 hisname.Form("hDLenBkgPt%dLS",i);
1171 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1172 fDLenHistLS[index]->Sumw2();
1173 hisname.Form("hSumd02BkgPt%dLS",i);
1174 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1175 fSumd02HistLS[index]->Sumw2();
1176 hisname.Form("hSigVertBkgPt%dLS",i);
1177 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1178 fSigVertHistLS[index]->Sumw2();
1179 hisname.Form("hPtMaxBkgPt%dLS",i);
1180 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1181 fPtMaxHistLS[index]->Sumw2();
1182 hisname.Form("hDCABkgPt%dLS",i);
1183 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1184 fDCAHistLS[index]->Sumw2();
1185
1186 indexLS++;
1187 hisname.Form("hLSPt%dLCntripsinglecut",i);
1188 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1189 fMassHistLS[indexLS]->Sumw2();
1190 hisname.Form("hLSPt%dTCntripsinglecut",i);
1191 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1192 fMassHistLSTC[indexLS]->Sumw2();
1193
1194 indexLS++;
1195 hisname.Form("hLSPt%dLCspc",i);
1196 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1197 fMassHistLS[indexLS]->Sumw2();
1198 hisname.Form("hLSPt%dTCspc",i);
1199 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1200 fMassHistLSTC[indexLS]->Sumw2();
1201 }
1202
1203 for(Int_t i=0; i<3*fNPtBins; i++){
1204 fOutput->Add(fCosPHistLS[i]);
1205 fOutput->Add(fDLenHistLS[i]);
1206 fOutput->Add(fSumd02HistLS[i]);
1207 fOutput->Add(fSigVertHistLS[i]);
1208 fOutput->Add(fPtMaxHistLS[i]);
1209 fOutput->Add(fDCAHistLS[i]);
1210 }
1211 for(Int_t i=0; i<5*fNPtBins; i++){
1212 fOutput->Add(fMassHistLS[i]);
1213 fOutput->Add(fMassHistLSTC[i]);
1214 }
1215}
1216
1217//________________________________________________________________________
1218void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){
1219 // Histos for impact paramter study
1220
1221 Int_t nmassbins=GetNBinsHistos();
1222 Int_t nbins[3]={nmassbins,200,fNImpParBins};
1223 Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
1224 Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
1225
1226 fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1227 "Mass vs. pt vs.imppar - All",
1228 3,nbins,xmin,xmax);
1229 fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1230 "Mass vs. pt vs.imppar - promptD",
1231 3,nbins,xmin,xmax);
1232 fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1233 "Mass vs. pt vs.imppar - DfromB",
1234 3,nbins,xmin,xmax);
1235 fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1236 "Mass vs. pt vs.true imppar -DfromB",
1237 3,nbins,xmin,xmax);
1238 fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1239 "Mass vs. pt vs.imppar - backgr.",
1240 3,nbins,xmin,xmax);
1241
1242 for(Int_t i=0; i<5;i++){
1243 fOutput->Add(fHistMassPtImpParTC[i]);
1244 }
1245}
1246
1247//________________________________________________________________________
1248void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1249{
1250 // Terminate analysis
1251 //
1252 if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1253
1254 fOutput = dynamic_cast<TList*> (GetOutputData(1));
1255 if (!fOutput) {
1256 printf("ERROR: fOutput not available\n");
1257 return;
1258 }
1259 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1260
1261 TString hisname;
1262 Int_t index=0;
1263
1264 for(Int_t i=0;i<fNPtBins;i++){
1265 index=GetHistoIndex(i);
1266
1267 hisname.Form("hMassPt%dTC",i);
1268 fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1269 }
1270
1271 TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1272 c1->cd();
1273 TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1274 hMassPt->SetLineColor(kBlue);
1275 hMassPt->SetXTitle("M[GeV/c^{2}]");
1276 hMassPt->Draw();
1277
1278 return;
1279}
1280//_________________________________________________________________________________________________
1281Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1282 //
1283 // checking whether the mother of the particles come from a charm or a bottom quark
1284 //
1285
1286 Int_t pdgGranma = 0;
1287 Int_t mother = 0;
1288 mother = mcPartCandidate->GetMother();
1289 Int_t istep = 0;
1290 Int_t abspdgGranma =0;
1291 Bool_t isFromB=kFALSE;
1292 Bool_t isQuarkFound=kFALSE;
1293 while (mother >0 ){
1294 istep++;
1295 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1296 if (mcGranma){
1297 pdgGranma = mcGranma->GetPdgCode();
1298 abspdgGranma = TMath::Abs(pdgGranma);
1299 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1300 isFromB=kTRUE;
1301 }
1302 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1303 mother = mcGranma->GetMother();
1304 }else{
1305 AliError("Failed casting the mother particle!");
1306 break;
1307 }
1308 }
1309
1310 if(isFromB) return 5;
1311 else return 4;
1312}
1313//_________________________________________________________________________________________________
1314Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const {
1315 // true impact parameter calculation
1316
1317 Double_t vtxTrue[3];
1318 mcHeader->GetVertex(vtxTrue);
1319 Double_t origD[3];
1320 partDp->XvYvZv(origD);
1321 Short_t charge=partDp->Charge();
1322 Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1323 for(Int_t iDau=0; iDau<3; iDau++){
1324 pXdauTrue[iDau]=0.;
1325 pYdauTrue[iDau]=0.;
1326 pZdauTrue[iDau]=0.;
1327 }
1328
1329 Int_t nDau=partDp->GetNDaughters();
1330 Int_t labelFirstDau = partDp->GetDaughter(0);
1331 if(nDau==3){
1332 for(Int_t iDau=0; iDau<3; iDau++){
1333 Int_t ind = labelFirstDau+iDau;
1334 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1335 if(!part){
1336 AliError("Daughter particle not found in MC array");
1337 return 99999.;
1338 }
1339 pXdauTrue[iDau]=part->Px();
1340 pYdauTrue[iDau]=part->Py();
1341 pZdauTrue[iDau]=part->Pz();
1342 }
1343 }else if(nDau==2){
1344 Int_t theDau=0;
1345 for(Int_t iDau=0; iDau<2; iDau++){
1346 Int_t ind = labelFirstDau+iDau;
1347 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1348 if(!part){
1349 AliError("Daughter particle not found in MC array");
1350 return 99999.;
1351 }
1352 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1353 if(pdgCode==211 || pdgCode==321){
1354 pXdauTrue[theDau]=part->Px();
1355 pYdauTrue[theDau]=part->Py();
1356 pZdauTrue[theDau]=part->Pz();
1357 ++theDau;
1358 }else{
1359 Int_t nDauRes=part->GetNDaughters();
1360 if(nDauRes==2){
1361 Int_t labelFirstDauRes = part->GetDaughter(0);
1362 for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1363 Int_t indDR = labelFirstDauRes+iDauRes;
1364 AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1365 if(!partDR){
1366 AliError("Daughter particle not found in MC array");
1367 return 99999.;
1368 }
1369
1370 Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1371 if(pdgCodeDR==211 || pdgCodeDR==321){
1372 pXdauTrue[theDau]=partDR->Px();
1373 pYdauTrue[theDau]=partDR->Py();
1374 pZdauTrue[theDau]=partDR->Pz();
1375 ++theDau;
1376 }
1377 }
1378 }
1379 }
1380 }
1381 if(theDau!=3){
1382 AliError("Wrong number of decay prongs");
1383 return 99999.;
1384 }
1385 }
1386
1387 Double_t d0dummy[3]={0.,0.,0.};
1388 AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1389 return aodDplusMC.ImpParXY();
1390
1391}