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